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PREFACE 


A paper  entitled  "Thermal  Response  Factors  for  Multi-layer  Structures 
of  Various  Heat  Conduction  Systems"  was  published  in  the  1969  Transactions 
of  the  American  Society  of  Heating,  Refrigerating  and  Air  Conditioning 
Engineers  (Vol.  75,  Part  1,  pp.  246-271).  A computer  program  mentioned 
in  that  paper  ailed  RESPTK  was  used  to  obtain  thermal  response  factors 
for  walls,  solid  objects  and  semi-infinite  walls  of  planes,  cylindrical 
and  spherical  systems.  In  order  to  respond  to  the  many  requests  for 
the  computer  program,  the  Fortran  listing,  the  input  instructions  of 
RESPTK  and  various  modifications  to  that  paper  have  been  included 
in  this  report. 

Contingent  upon  the  responses  given  to  this  report  a formal  NBS 
publication  may  be  issued  in  the  future  as  a Building  Science  Series 


paper . 


ABSTRACT 


The  thermal  response  factor  method  for  calculating  transient  heat 
conduction  through  multi-layer  slabs  is  generalized  to  include  the 
solutions  for  many  other  important  engineering  heat  transfer  problems. 
Response  factor  formulas  for  multi-layer  structures  of  cylindrical 
and  spherical  objects  (hollow  as  well  as  solid),  plane  and  curved 
surface  walls  adjacent  to  an  infinitely  thick  heat  conduction  medium, 
such  as  the  ground,  and  to  plane  slabs,  are  presented  in  this  paper. 
Numerical  evaluation  of  these  formulas  is  carried  out  for  selected 
multi-layer  structures  and  results  are  tabulated. 

Also  included  in  this  report  are  Fortran  listings  of  the  response 
factor  calculation  programs  and  sample  usages  of  the  computer  programs 
for  evaluating  heat  conduction  through  building  walls. 

Key  Words:  Thermal  response  factors,  multi-layer  structures, 

transient  heat  conduction,  cylinder  and  sphere 


Nomenclature 


Unless  otherwise  specified,  the  following  symbols  are  used  through- 
out this  paper.  Since  the  units  attached  to  symbols  represent  the 
English  system  (still  most  popular  among  heat  and  air  conditioning 
engineers  in  the  United  States),  a conversion  table  for  standard  metric 
units  is  also  provided  at  the  end  of  this  section. 


A , B , C , D 


A , B , C , D 

V*  V V V 


hI*  ho 


J 

X 


V 


m 


n 


Elements  of  overall  temperature  flux  matrix 
Elements  of  individual  layer  temperature-flux 
matrix 

-1  -2 

Heat  flux,  Btu  hr  ft 
Laplace  transform  of  heat  flux 
Exterior  wall  surface  heat  conductance 
[Btu  ft"2  F"1  hr"1] 

-2  -1  -1 

Irradiated  heat  flux,  [Btu  ft  hr  F ] 


Complex  number  notation 


= /T 


Thermal  conductivity  of  the  vth  layer 
[Btu  hr"1  F"1  ft"1] 

Thickness  of  the  vth  layer,  fft.] 

Curvature  index:  m = 0,  plane;  m = 1,  cylin- 

der; m = 2,  sphere 
A large  number 

Total  number  of  layers  to  be  considered 
Laplace  transform  parameter  (p  is  treated  as  a 

complex  variable  for  the  inversion  integral) 

-1  -2 

Heat  flux,  fBtu  hr  ft  ] 


i 


q 

ii 

e roj 

r , r , 

v*  v+1 

Radii  of  the  bounding  surfaces  of  the  vth 

R 

V 

layer , [ft] 

2 

Thermal  resistance  of  the  vth  layer,  [ft  F 
hr  Btu  ^] 

R 

General  response  function  defined  in  the  text 

T 

Temperature,  [F] 

T 

o 

Initial  temperature  at  t = 0,  fF] 

T , T . 
v v+1 

Boundary  temperatures  of  the  vth  layer,  fF] 

V , V - 

v v+1 

Temperature  departure  of  the  vth  layer,  (T  - 

V 

Tq)  and  (T^^-T^),  respectively,  TF] 

V , V , 

v v+1 

Laplace  transforms  of  the  temperature  depar- 

Xf» and  Z_^ 

tures,  V and  V , 
s v+1 

-2  -1  -1 

Response  factors,  [Btu  ft  F hr  ] 

Xt,  Y1#  and  Z. 

-2  -1  -1 

Modified  response  factors  [Btu  ft  F hr  ] 

% 

2 -1 

Thermal  diffusivity  of  vth  layer,  [ft  hr  ] 

«P 

Characteristic  function  of  pulses 

6 

Time  increment,  fhr] 

A 

V 

Determinant  of  the  matrix  for  A , B , C , D 

V V V V 

r 

Determinant  of  the  matrix  of  A,  B,  C,  D 

t . 

i 

Time  coordinate  t^  = t - i§,  [hr] 

T 

Time  index  t = t6 

Q 

p* 

3^  = -p^»  roots  for  residue  evaluation 
k = 1,  2,  3,.... 

Defined  in  the  text,  eq.  (26) 
ii 

Subscripts 


v = layer  boundaries 

■Pt  = roots  for  the  residue  evaluation 

i = response  factor  series  in  relation  to  time  series 
t = discrete  time 


Unit  Conversion 


To  convert  from 

Multiply  by 

To  Obtain 

Q = Btu  hr"1  ft"^ 

3 . 152481E  + 00 

W m"2 

h = Btu  hr-1  ft"2  °F” 1 

5.6783E  + 00 

W m" 2 K"1 

\ - Btu,  hr"1  ft"1  0F-1 

0. 14423E  - 02 

W m“ 1 K"1 

c * Btu  lb-1  F"1 

4.187E  + 03 

J Kg"1  °K" 1 

p = lb  f t”3 

1.602E  + 01 

Kg  m"  3 

a = ft2  hr"1 

2.581E  - 05 

m2  s" 1 

l = ft 

3.048E  - 01 

m 

iii 


1.  Introduction 


A recent  advance  in  computer  application  for  hour  by  hour  building 
heat  transfer  calculations  has  made  it  possible  to  improve  the  Response 
Factor  technique  in  transient  heat  conduction  analysis.  This  improved 
response  factor  method  permits  an  accurate  evaluation  of  transient  (non- 
steady and/or  aperiodic)  heat  conduction  through  multi-layer  walls  and 
roofs,  which  has  heretofore  been  extremely  difficult. 

It  should  be  stated  that  an  existing  procedure  commonly  known  in 

1-2/ 

the  U.  S.  as  the  Mackey  and  Wright — — solution  for  evaluating  the 
building  heat  transfer  has  been  based  upon  the  assumption  that  the 
building  walls  and  roofs  experience  steady  periodic  temperature  cycles 
on  a diurnal  basis.  Their  solution  obtained  using  this  assumption  is  in- 
adequate for  the  accurate  evaluation  of  actual  hour-by-hour  heat  gain 
or  loss  of  buildings. 

Another  well-known  approach  solving  the  transient  heat  transfer  is 
finite  difference  approximations  to  the  heat  conduction  equation.  Al- 
though the  computational  procedures  involved  in  this  latter  technique 
are  less  complicated  than  analytical  procedures,  extremely  small  grid 
sizes  are  required  for  finite  difference  time  and  space  coordinates  if 
computational  stability  is  to  be  retained  in  calculating  transient 
heat  flow  for  a multi-layer  heat  flow  problem. 


The  Response  Factor  Method  has  been  treated  previously  by  several 

30/ 

authors—-— 7 . This  method  basically  utilizes  the  superposition  principle 
in  such  a manner  that  the  overall  thermal  response  of  the  building 
structure  at  a selected  time  is  the  sum  of  the  responses  caused  by  many 
individual  temperature  pulses  during  preceding  significant  times.  Thus, 
by  simulating  the  transient  boundary  temperatures  by  a train  of  pulses, 
and  by  summing  up  the  heat  flux  caused  by  each  pulse,  the  total  heat 
flux  at  a given  time  can  be  derived.  The  calculation  of  the  thermal 
response  of  multi-layer  walls  and  roofs  due  to  each  individual  tem- 
perature pulse  has  in  the  past  been  simulated  by  the  concept  of  a finite 
number  of  lumped-resistances-and-capacitances  by  equating  the  heat  flow 

path  to  an  electrical  circuit  analog*  A significant  contribution  has 

9/ 

been  made  recently  by  Mitalas  and  Arseneault—  , who  improved  considerably 

the  accuracy  of  the  calculation  by  avoiding  the  lumped-resistance-and- 

capacitance  concept.  Mitalas  and  Arseneault  were  able  to  solve  the 

differential  equation  of  heat  conduction  for  the  multilayer  system  by 

employing  a matrix  equation  of  Laplace  transforms—^.  Although  the 

matrix  equation  of  the  composite  wall  has  been  employed  effectively  in 

1 2 6 7-8/ 

the  past,  previous  efforts — 2 — *— ■*■  have  used  Fourier  series  simulation 

of  the  boundary  temperature  functions.  The  previous  difficulty  of 

applying  the  matrix  method  for  a periodic  heat  transfer  problem  was 

1 

\ 

primarily  due  to  the  complexity  of  evaluating  the  inversion  integrals 
of  the  Laplace  transforms.  Using  a high-speed  digital  computer,  Mitalas 


2 


and  Arsenault  were  able  numerically  to  invert  the  Laplace  transform 
matrix  for  the  multi-layer  heat  conduction  equation  when  a periodic  or 
transient  boundary  temperature  function  is  simulated  by  a train  of 
triangular  pulses. 

In  this  paper,  the  method  employed  by  Mitalas  and  Arseneault  is 
extended  to  cover  walls  and  roofs  with  cylindrical  and  spherical  curva- 
tures. The  response  factors  for  the  cylindrical  and  spherical  walls 
may  also  be  useful  in  analyzing  the  transient  heat  flow  through  pipes, 
underground  shelters  and  structures,  storage  tanks,  and  tunnels.  Sample 
calculations  for  a typical  brick  wall  were  performed  using  a computer 
program  to  carry  out  the  mathematical  procedures  outlined  in  this  paper. 
Results  of  the  calculations  for  plane,  cylindrical  and  spherical  walls 
are  compared  with  those  obtained  by  an  exact  analytical  method  for  a 
steady  periodic  boundary  temperature  profile. 

Also  presented  in  this  paper  are  formulas  for  evaluating  heat  flux 
in  semi-infinite  systems,  interfacial  temperatures  and  heat  fluxes  of 
the  multi-layer  transient  system,  and  the  analyses  for  non-linear  bound- 
ary heat  transfer  problems. 
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2 . Heat  Conduction  Through  a Homogeneous  Layer 


The  heat  conduction  equations  for  one -dimensional  heat  flow  in  a 
homogeneous  layer  of  a multilayer  system  are  first  analyzed.  Assume 
that  this  particular  layer  has  thermal  diffusivity  thermal  con- 

ductivity and  at  time  t has  boundary  temperatures  T^(t)  at  the 

surface  r - r . and  T . (t)  at  the  surface  r - r , . Also  assume  that 
V*  V+l  v+l 

the  temperature  of  the  layer  at  time  t = 0 was  constant  at  T . The 

o 

differential  equation  and  boundary  conditions  describing  the  conditions 
stated  above  are  then 


S2T  m ST  _ ST 
Sr2"  + r Sr  St 


for  m = 0,  1,  or  2 


T = T at  r = r 

v v 


T = T , at  r = r , 

v+l  v+l 


T = T for  all  r at  t = 0 
o 


for  t > 0 


(1) 


Applying  the  Laplace  transform  to  the  above  relations,  it  is  possible 
to  write 

V = V at  r = r 

v v 


d fV 

dr® 


m dV 

r dr 


q..v 


(2) 


V = V . at  r = r , 

v+l  v+l 


where 


V = f (T  - T0)e-Pfcdt 


qv  = rr 

V Qv 


(3) 

(A) 


and  p is  the  Laplace  transform  operator. 
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A general  solution  of  these  Laplace  transform  equations  for  heat  conduc- 


tion may  be  written  in  matrix  form, 


d'OGM 

V v+1 


(5) 


or  by  rearranging. 


(6) 


v 

dV 

where  f and  f , are  the  Laplace  transforms  of  — , heat  flux  at 
v v+1  vdr  * 

r = r , and  respectively.  Specific  expressions  for  each  element 

of  the  matrix  in  (5)  for  the  cases  of  m = 0,  1,  and  2 are  shown  in 
Tables  1,  2 and  3. 

By  using  the  expressions  in  Tables  1,  2,  and  3,  it  can  be  shown 
that  the  determinant  of  the  matrix  in  (5)  is 


r 

V 


A 

CV 

v 


(7) 


3.  Multi-layer  Heat  Conduction 

The  solutions  obtained  for  the  single  layer  (the  vth  layer)  bounded 
by  r = r^  and  r = r^+1  are  valid  for  each  of  the  other  layers  of  a multi- 
layer slab,  so  that  one  may  write  for  each  layer: 


1st 

layer : 

G;)  = 

(Al 

V/M 

V ij 

2nd 

layer : 

(+ 

VCj 

D8Af3/ 

(8) 

(n-1) st 

layer : 

(Vn-1) 

_ fAn 

- 1 ®n- 1 ^n 

V fn-l/ 

-1  Dn-1/Afn/ 

5 


This  is  predicted  upon  the  assumption  that  there  is  perfect  thermal 
contact  at  the  interface  of  the  layers  of  the  multi-layer  slab  giving 
continuity  of  temperature  and  heat  flux.  Combining  the  above  matrix 
equations  give 

'n 

where 


(y  --  ft  ®YV) 

\fn  / \C  D/\,.  J 


(9) 


'A 

B\ 

- 

V) 

f A2 

V)  ( 

\C 

VC, 

Di  / 

\CQ 

A i B 

Cnft  D n-]J 

n-1  n-1 


(10) 


If  the  vth  layer  of  the  multi-layer  slab  has  a negligibly  small  thermal 
mass,  (e.g.,  is  a fully  enclosed  air  space),  the  matrix  elements  for  that 
layer  are 


K = i 

Bv  = Rv 


(11) 


cv  = 0 

Bv  - rv 

where  R is  the  thermal  resistance  of  the  layer, 
v 3 

Applying  matrix  algebra,  the  determinant  of  the  overall  matrix 
in  (10)  can  be  shown  to  be 


A B 
C D 


= r = r,  . ra  • r. 


r - 

..  _ !^vj 


m 


(12) 


From  (9),  the  Laplace  transform  of  the  heat  flux  matrix  relation  is 


\ d _r  _ 

1 j = {B  B'l/Xl  ) 

V,  AV 

JL  _ A 
B ~B 


(13) 


The  heat  flux  at  each  surface  can  be  evaluated  by  applying  the  inversion 
theorem  of  the  Laplace  transform  to  equation  (13). 
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4.  Superposition  Principle  and  the  Inversion  of 
Laplace  Transforms 

The  inversion  of  (13)  can  be  approximated  easily  by  applying  the 
superposition  principle  where  the  slab  temperature  T is  represented  by 

5 a linear  sum  of  functions  V^(i  = 1,  2,  ....)  such  that 

CO 

T - T0  = r Vi(ti)  (14) 

i=0 

Furthermore,  the  boundary  temperature  functions  Tx  and  Tn  at  r = rx 
and  r = rn  are  assumed  to  be  represented  by  a series  of  pulse  functions 
n such  that 


Vi 


Ti 


^ Vi,i<p(ti) 

i=0 


V = T„  - T = 
n n xo 


Vn,i  <p(ti> 


i=0 


15  In  the  above  equation,  V^j  and  are  pulse  heights  at  time  t = i5 

for  the  boundary  surfaces,  where  6 is  the  discrete  time  interval  of  the 
pulses.  The  pulse  function  cp(t^)  is  defined  only  for  0 < t^  < m'6} 
where  m'  is  the  width  of  the  pulse  at  the  time  base.  The  simplest  pulse 
most  commonly  used  is  the  rectangular  pulse  of  width  5,  (or  m'  =1), 

20  such  as  shown  in  Fig.  1,  and  it  can  be  described  by  the  following 
pulse  function 

cp ( 1 1 ) =0  ti  < 0 

= 1 0 < ^ 6 
=0  ti  > 6 
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Fig.  2,  Trapezoidal  pulses 


Although  the  rectangular  pulse  simulation  of  the  boundary  temperatures 
is  very  simple,  the  approximation  of  a complex  profile  by  a finite 
number  of  rectangular  pulses  inevitably  causes  loss  of  accuracy  unless 
the  time  increment  6 is  chosen  extremely  small.  A considerable  gain  in 
the  accuracy,  however,  can  be  restored  if  the  boundary  temperatures  are 
simulated  by  trapezoidal  pulses,  such  as  shown  in  Fig.  2.  It  can  be 
proven  also  that  two  overlapping  triangular  pulses  (dotted  line)  of 
base  width  of  26  have  identical  thermal  response  to  that  created  by 
the  trapezoidal  pulse  of  width  6,  which  is  shared  by  the  two  triangular 
pulses  (see  Fig.  2).  The  triangular  pulse  of  m'  =2  is,  however,  better 
suited  for  this  analysis  than  the  trapezoidal  pulse,  since  it  represents 
each  pulse  by  a single  pulse  instead  of  two.  The  pulse  function  for  a 
triangular  pulse  of  base  26  is 


o 

II 

✓ — - 
•r-i 
■U 

for 

t . 

l 

0 

= t./6 

i 

for  0 < 

t. 

i 

6 

(16) 

= 2-  t./6 
i 

for  6 < 

t . 
l 

26 

= 0 

for 

t. 

> 

26 

Substituting  (14)  and  (16)  into  the  original  differential  equation  it 

is  found  that  the  solutions  obtained  for  V are  also  valid  for 

V.  (i  = 1,  2,  ....  co)  provided  that  the  new  time  coordinate  t.  used 
i i 

for  V_j.  is  related  to  the  original  time  coordinate  t by 

t.  = t - i6 

l 


i 

! 
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The  Laplace  transform  flux  relation  for  V,  . and  V . , similar  to 

l,i  n,i» 

equation  (13)  is  then  written  as 


d _ r 


(S:i)  ■ * (?  B)(^>t-i6) 

1 - y^vn,t-i57 


(17) 


B B 

where  cp  is  the  Laplace  transform  of  the  pulse  function  cp,  or 

¥ = for  0 < fci  * 6 


~ 6p^"  “ 2e~P&)  for  6 < tt  ^ 26 

= 5^5-  (1  - e“P6)2,  for  26  < t± 


(18) 


for  the  triangular  pulse  function. 

The  inversion  of  the  Laplace  transform  can  be  accomplished  by 
applying  the  residue  theorem  to  the  inversion  integral,  details  of  which 
i u 10/ 

are  given  elsewhere — . 

The  inversion  of  flux  equation  (17)  essentially  involves  the 

analysis  of  the  following  general  formula 

_ R 

fi  = cp  B (19) 

where  R represents  D,  T,  1,  or  A in  equation  (17).  The  inversion  of 


(19)  yields 


d_  rp2cpRepti-i  r“  rcpRe 
'l  * Urn  dp  . ! + 'S  ^ 


Pfci 


1 


p-*0 


ft=l 


L 8b 

dp 


JP=-8ft 


(20) 


where  3 (ft  = 1,  2,  ....)  are  the  real  roots  of  the  equation 
ft 

B(p  = -8^)  = 0.  (21 

In  order  to  evaluate  equation  (20),  values  of  A,  B,  D,  ^ at 

p = 0 and  p = -3 . , (ft  =0,  1,  2,  ....)  are  needed, 
ft 
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The  derivatives  of  matrix  elements  A , B , C , and  D are  provided 

V V V*  V 

in  tables  4,  5,  and  6 to  assist  in  the  calculation  of  the  derivatives  of 
A,  B,  C,  and  D by  the  following  relationship: 


d_  /A 

dp  Vc 


BN 

DV  ~ 


+ 


+ 


dAx 

dq 

/dp 

dp  )($ 

^dq 

dDx 

dp 

dp 

/Ai 

B \ dA^ 

bl 

v p 

vq 

XACz 

dp 

dB3 

dp  \ 

)• 

£V 

dp 


/ Ai  Bx  / Ag  B2\ 

vCi  DtAOs 


/ An- 1 
N^n-1 


dAn-l 
Z'  dP 

' dCn-l 
dp 


Bn-1 

Dn-Iy 


Bn-1\ 


dBn-1 

dp 

dPr-/ 

dp 


(22) 


Applying  the  theory  of  limits. 


it  can  also  be  shown  that 


*“  V,  _ A V, 
P-°  \Cv  I v "Vo  1 J 


V 


.1 

/ 

\ 

0 


^v+l\ 


v+lY3 


A 

) 


t 

V 


for  m = 0 


for  m = 1 (23) 


for  m = 2 
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It  is  extremely  interesting  to  note  that  the  multiplication  of  these 
successive  matrices  for  the  multi-component  slab  would  yield 


, . 1 
Urn  /A  B\  = /I  T7y 

P-^  \C  D J V0  r J 


where  U'  = overall  steady-state  heat  conductance  from  the  surface  1 to 
n + 1. 

The  first  term  of  the  right-hand  side  of  equation  (20)  can  be 
reduced  further,  for  the  case  of  the  triangular  pulse  function,  to  the 
following  relation 


1 im 

p-*  0 


d_  r p2cpRePti-| 
dP  LB ' 


U'rdR  + 
6 ( dp 

U'rdR 
6 !_dp 


R ^ B-, 
R - R dp 


B JP=0 , 


R — 1 


dp 


P=0 


= 0 


J 


0 < t^  s'  § 

6 < t.  <;  26 
26  < tt 


(24) 


Table  7 is  provided  for  the  evaluation  of  the  limits  of  each  of  the 

dR 

matrix  elements  and  their  derivatives,  such  as  R and  — , as  p approaches 


zero. 


Letting  p = -p  , then 

n 


= J 


1 3^ 


(25) 


whereby  the  complex  functions  of  Tables  1,  2,  3,  4,  5,  and  6 and  the 
derivatives  occurring  in  series  on  the  right-hand  side  of  equation  (20) 
can  be  represented  as  real  functions  as  shown  in  Tables  8,  9,  and  10. 
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The  functions  indicated  in  Tables  8,  9 and  10  are  to  be  evalua- 
ted at  each  of  the  negative  real  roots  -{3^  (ft  = 1,  2,  3,  ) of  the 

equation  B(p)  = 0,  for  all  non-negligible  terms  of  equation  (20).  The 
magnitude  of  the  terms,  however,  decreases  quite  rapidly  with  increase 

in  ft,  particularly  when  t.  is  large  or  when  a particular  component  has 

1. 


a large 


/o \ 


value . For  the  triangular  pulse  function,  the  series  of 


equation  (20)  can  be  evaluated  by  the  following  relations 


Y rcpRepti  - 

cfB  ~ 


ft=l 


dp 


p=-bft  $=i 


' *(|3ft)  e 


-n 


for  t^  < 6 


Y(8ft)(l  - 2eX)  e"2Q 


ft=l 


for  6 < t^  ^ 26 

(26) 


T (P*)(l  - eX)s  e 


■in 


ft=i 


for  t£  = i6  > 26 


where 


1 r? 


- W l!5  . 


dp-  p=- 


Cl  — 8^6 


where  R may  be  any  one  of  A,  F,  1,  or  D of  Equation  (17) 
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By  combining  (24)  and  (26),  generalized  response  factors 

(i  = 0,  1,  2,  ....  a>)  may  be  derived  in  terms  of  R and  its 

derivative  as  follows: 
dp 


Xn  = m 


^dR 
+ dp 
i_B_ip=o  L'B6 


R? 

B^6  -'p=0 


+ (for  1 

^=1 


= 1) 


X 


1 


pdR  dB 

"j  dp  - R Bp  1 
LB6  B2  6 P=0 


+ ^ Y(3ft)(i-2e0^6)e-20^6 

fr=l 


(for  i = 2) 


(27) 


V Y(Pft)  (l-e^6)2  e_i^6 

Ui 

*=l 


(for  i=3,  4,...®) 


Using  these  notations,  the  inversion  of  heat  flux  relation  (17)  may  be 
expressed  generally  as 


/Fl,i\  = /X  i -rYiYVl(t-i6)' 
^Fn,i;  KYi  "zi  Avn(t-i6); 


(28) 


where  X^,  Y^,  and  are  response  factors  and  correspond  to  X^  of 

equation  (27),  (28)  and  (29)  when  R is  replaced  by  D,  1,  and  A respectively. 

By  denoting  the  time  coordinate  t by  increments  of  6,  say  t = t6, 

V1(t  - i6)  may  be  expressed  simply  by  V1  ..  Using  the  subscripted 
JL  1 y T“  1 

temperature  notation,  equation  (28)  can  be  used  to  express  the  original 
heat  conduction  system  as  follows 


(F1,A 

VFn,r2 


^ (Xi  -rYY(/Vi}T-i' 
L*  Ui  -Zi  /\Vn,T-i^ 
i=0 


(29) 


This  relation  is  called  the  convolution  equation  of  the  heat  fluxes. 
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In  equation  (29),  X^  Y±,  and  Z.  are  called  response  factors.  A 
close  examination  of  (27),  (28)  and  (29)  reveals  the  following  interesting 
facts: 


1.  Response  factors  and  Z±  tend  to  decrease  with  a common 

R 8 

ratio  e Pl  for  large  values  of  i or 


Xi+1  — Yi~H  = zi+l  = if  i > N (30) 

Xi  Yi  Zi 

and  N is  a large  number.  For  a conventional  building  wall, 

N ss  15. 

2.  To  be  compatible  with  the  steady  state  heat  flow  condition  when 

V.  and  V are  constant,  it  is  necessary  that 
1 n 


CD 

00 

CO 

'T"’  ,, 

v 1 

L Xi 

ss 

Yi 

= 

) zi 

L — i 

— 1 

i~0 

i=0 

i=0 

(31) 


where  U is  the  overall  heat  transfer  coefficient  based  upon 

r = rn+r 

5.  Sample  Calculations 

A digital  computer  program  called  RESPTK  (refer  to  Appendix)  has 
been  developed  at  the  National  Bureau  of  Standards  for  calculating  the 
response  factors  formulated  in  the  previous  sections.  Sample  walls  with 
properties  as  shown  in  Fig.  3 and  Table  11  were  analyzed  by  this  program 
for  cases  m = 0,  1,  and  2 (for  plane  wall  (PW),  cylindrical  wall  (CW), 
and  spherical  wall  (SW),  respectively).  The  sample  wall  consists  of 
two  solid  mass  layers  bounded  by  two  air  film  layers. 
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Table  12  shows  the  residues  of  _Jp,  and  _J>  at  p = 0 for  0 < t ^ 6 

B B B 

and  for  6 < t ^26.  The  residue  of  these  functions  becomes  zero  for 
t > 26. 

Table  13  gives  3^,  the  roots  of  B(p)  = 0 along  the  negative  real 
axis.  The  response  factors  calculated  by  formula  (20)  for  R = D,  1 and 
A are  indicated  in  Table  14  as  X^,  Y^,  and  Z^,  respectively.  Also 
indicated  at  the  end  of  Table  14  are  the  common  ratios  from  (30) 
attained  by  successive  values  of  each  of  the  response  factors  when 
i > 14.  Each  of  the  response  factors  corresponds  to  the  value  evaluated 
at  t = i5.  As  seen  from  Table  14,  the  response  factors  for  plane,  cylin 
drical,  and  spherical  walls  are  very  similar  to  each  other  in  this 
particular  wall.  This  is  due  to  the  fact  that  the  curved  walls  used 
in  the  sample  calculations  had  an  innermost  radius  of  5 ft.  and  total 
wall  thickness  of  2/3  ft.,  which  can  be  very  closely  simulated  by  the 
plane  wall  heat  transfer.  Using  the  response  factors  the  heat  flux 
values  at  r = r^  and  r = r^  were  also  calculated  for  a periodic  tem- 
perature profile,  results  of  which  are  shown  in  Table  15,  16,  and  17 
corresponding  respectively  to  the  plane,  cylindrical  and  spherical  walls 
Although  the  application  of  the  response  factor  calculation  is  not 
limited  to  the  periodic  heat  flow  problem,  the  periodic  heat  flow  pro- 
blem was  chosen  for  the  sample  calculation  because  exact  solutions  for 
the  periodic  heat  flow  problem  can  be  used  to  check  the  accuracy  of  the 
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response  factor  method.  (The  response  factor  calculation  is,  in  a rigo- 
rous sense,  an  approximate  solution,  where  the  boundary  temperature 
profiles  are  approximated  by  a train  of  trapezoidal  pulses).  The  exact 
solutions  for  the  heat  conduction  equation  under  periodic  boundary 
temperature  conditions  are  obtained  by  setting  in  the  original  differen- 
tial equation 

00 

Tx  = T0  + V V1}.  eju)it 

i=l 


T = 


T + 
Ao 


I Vn’ 

i=l 


>ifc 


(34) 


qv,i  = /fE  <1+J)>  »iafa”d?  = 24 

The  heat  flux  relations  in  terms  of  complex  variables  are  treated  in 
reference  [10].  Another  computer  program  called  ETD  2 was  developed 
during  the  course  of  this  study  to  perform  the  complex  algebra  calcula- 
tion for  the  periodic  heat  flow  problem.  The  results  of  the  exact 
solutions  are  given  in  Tables  15,  16,  and  17. 

Agreement  between  the  exact  solutions  and  the  solutions  obtained 
by  the  response  factor  calculations  shown  in  Tables  15,  16,  and  17  is 
very  good.  To  obtain  this  degree  of  agreement,  the  response  factors  had 
to  be  calculated  up  to  i = 72 . The  compilation  and  computation  of  heat 
flux  by  response  factors  for  all  three  walls  using  a UNIVAC  1108  took 
34  seconds,  while  for  a periodic  heat  transfer  solution  by  complex  alge- 
bra, the  time  was  16  seconds.  The  response  factor  calculation  involves 
a lengthy  iterative  process  in  searching  for  the  roots  of  B(p)  = 0. 
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The  plane  wall  response  factors  treated  in  this  paper  have  also 
been  calculated  by  D.  G.  Stephenson  of  the  National  Research  Council  of 


For  the  calculation  of  heating  and  cooling  loads  for  buildings,  it 
may  become  necessary  to  determine  the  heat  storage  effect  of  interior 
furnishings,  partitions,  floors,  ceilings,  etc.  This  heat  exchange 
problem  may  also  be  treated  by  the  response  factor  method  if  these 
materials  can  be  represented  by  simple  geometric  shapes  such  as  solid 
plates,  cylinders  or  spheres. 

c)T 

The  boundary  conditions  for  these  cases  are  — = 0 at  r = 0 for 


all  t > 0 for  the  cylindrical  and  spherical  cases  (m  £ 0).  For  the 
case  of  cylindrical  and  spherical  objects,  the  Laplace  transform  heat 
flux  relation  at  the  outside  surface  of  the  innermost  core  r = r1  is 


Canada- 


11/ 


His  results  agree  very  well  with  those  obtained  in  this 


paper . 


6.  Heating  and  Cooling  of  Plates,  Cylinders  and  Spheres 


expressed  in  terms  of  that  core's  thermal  properties  and 


fx  = G'Vi 


(33) 


where 


(34) 


(35) 
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Combining  this  relationship  with  the  rest  of  the  outer  multi-layer  system 
as  before 


Uny 


/D 

* 

i 

B 


. lAV 

B 


(36) 


and  noting  that 


d - r - 

fi  = B * B Vn 


G’Vj.  , 


(37) 


Vi 


vIT-  G' 

B 


(38) 


then  fn=^L-Av  (39) 

B B n 

. AG*  - C 77 
D - BG'  n 


The  inversion  of  this  heat  flux  relation  is  readily  obtained  by  the  resi- 
due theorem  similar  to  equation  (20).  Table  18  shows  specific  expressions 
of  G'  dG'/dp  for  m = 1 and  2 (or  for  the  cylinder  and  sphere. 

For  a plane  shaped  object  heated  or  cooled  at  both  surfaces  in  a 
space,  the  response  factor  representations  of  heat  exchange  between  the 
objects  and  the  air  in  space  at  time  t = t5  is 


CD 


qT  = ) (Xi  + Z.  - 2Y^  Tt_ i 


(40) 


i=0 

where  X_^,  Y^,  (i  = 0,  1,  2...)  are  response  factors  of  the  slab  as 
defined  in  (30)  including  surface  heat  transfer  coefficients,  and  T__  ^ 
represents  space  air  temperature  at  time  (t  - i)5. 
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7.  Semi- inf inite  System 

In  many  cases  the  transient  heat  conduction  characteristics  of 
semi-infinite  and  composite  systems  are  needed.  Problems  of  heat  conduc- 
tion to  the  paved  earth  surface,  to  the  basement  floor  and  to  underground 
pipes  are  good  examples  for  the  semi-infinite  system. 

Assume  the  nth  layer  of  the  previous  system  (equations  8,  9,  and 
10)  to  be  infinitely  thick,  its  thermal  conductivity  and  diffusivity 
been  and  . For  the  infinitely  thick  nth  layer,  boundary  conditions 
can  be  written  as  follows 


T = Tn  (t)  at  r = rn 
T = Tq  for  all  t at  rn+^  -*  oo 


(41) 


The  general  solution  in  the  Laplace  transform  domain  is 


fn  = GVn 


(42) 


where 


G = ^n<3n  for  m = 0 

= x n 

LSq,0"J^»r  m " 

~ 1 + V ~ ,for  m = 

l_  4n  Ln  J 


1*  =p= 


(43) 


Combining  relation  (37)  with  (9),  the  Laplace  transform  of  the  heat 
flux  equation  at  r = r,  (or  at  the  surface)  can  be  written  as 


f,= 


C_4_DG> 
\A  + B Gj 


Vi 


(44) 
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The  inversion  of  (44)  cannot  be  obtained  by  the  residue  theorem  as  in 
the  previous  cases  because  of  a branch  point — ^ at  p = 0.  The  branch 
point  integration  described  in  reference  flO]  may  be  performed  to  yield 
the  following  response  factor  relations  for  the  surface  heat  flux  over 
a multi-layer  semi-infinite  system. 

00 

= / ^i  Tl,T-i 

i=0 


where  Zq  = 

Zl  = 0g  “ 20j 

Zi  = 0i  - 20^„^  + 0j__2  for  i s 3 

®i  - r1V^i_  r (c  + dg)  szi  dp'>  dt 

o ' 2^ i J d and  c2  ^ ^ 


(45) 


and  cx  and  c2  are  paths  of  the  line  integral  where  p is  defined  by 
p = re'*"17  and  p = re  1TT  respectively  for  r from  zero  to  infinity.  Al- 
though the  method  described  above  is  for  a rigorous  and  generalized 
evaluation  of  response  factors,  an  approximate  solution  to  the  problem 
may  be  obtained  as  described  below. 

For  an  approximate  method,  the  response  factor  calculations  will 
be  performed  for  the  surface  layer  region  (which  may  or  may  not  be  of 
a multi-layer  system)  and  for  the  semi-infinite  region  separately. 

The  response  factors  for  the  former  are  denoted  herein  by  X^,  Y ^ and 
Z±  (i  - 0,  1,  2,  ....)  and  for  the  latter  by  Zl^  (i  = 0,  1,  ). 
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The  surface  temperature  of  the  entire  region  and  the  interfacial  tem- 
perature (at  r = rn ) between  the  surface  layer  and  the  semi-infinite 
regions  at  time  r are  denoted  by  ^ and  ^ respectively.  At  the 
interface,  the  following  heat  transfer  relations  can  be  established: 


n , t \ 7 1 fj, 

/ Li  rn,T-i 


/ 

i=0 


L YiTl,r-i  ~ ^ ZiTn,T-i 

i=0  i=0 


(46) 


Eliminating  ^ , from  the  above  equation,  and  applying  it  to  the 
heat  flux  at  the  surface  (r  = r1). 


1»t 


00 

v1  ■ 

L 

i=0 


Zi  T1 , T-i 


where  (i  = 0,1,... ®)  are  the  response  factors  for  the  overall  system 
(multi-layer  surface  region  plus  the  semi-infinite  region),  expressed  as 


Yj5 

Zi  + Zi  ' 


(47) 


Table  19  shows  mathematical  formulas  for  obtaining  Zl^  for  m = 0, 

1,  and  2.  Also  shown  in  the  Appendix  are  system  schematics  and  response 
factors  of  various  heat  conduction  systems  treated  in  this  paper. 
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8.  Interfacial  Temperatures 


The  response  factor  method  can  also  be  used  to  calculate  the  tem- 
peratures at  the  interfaces  of  a multi-layer  wall  system  when  surface 
temperatures  at  the  outer  surfaces  are  prescribed. 

For  the  multi-layer  system  of  Fig.  3,  assume  that  the  interfacial 
temperature  at  r = r3  is  to  be  calculated.  The  Laplace  transforms  of 
the  matrix  relations  for  the  flux  and  temperature  for  the  two  subsystems, 
one  for  the  layers  1 and  2,  and  the  other  for  the  layers  3 and  4,  are 
written  as  follows: 


1 


1 


By  denoting 


b(i)  \ 
Dd  )J 


and  by  knowing  that 


T\  1 


and 


B B 
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it  can  be  shown  that 


(48) 


The  inverse  transform  can  be  carried  out  to  yield  the  following 
relation,  if  the  triangular  pulse  simulations  are  used  for  and  V . 


3>t  ' £ (ai  VI,r-i  + bi  vo,T-i> 

i=0 


(49) 


where,  for  example,  b^  can  be  evaluated  by 


b. 

1 


d b(2)  - ip6H  V I b(2)  cpe  /CAX 

= Inn  — — — — cpe  + / — L (50) 

dp  _ B6  ^ J L L a.  dB  Jp=p . 


p-0 


fi=l  P 6 dp 


when  p is  the  ftth  negative  real  root  of  B(p)  = 0. 

7? 

9.  Application  of  Response  Factors  Calculation  to 
hlcf,  - Lirvw  Boundary  Problems 

In  many  heat  conduction  problems,  non-linear  heat  transfer  relations 
occur  at  boundary  surfaces.  Two  such  cases  of  major  importance  are 
treated  in  this  section  as  illustrative  examples  of  the  response  factors 
technique. 

Case  1:  Stefan-Boltzmann  type  radiation  heat  exchange  at  one  of  the 

surfaces: 
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This  situation  is  typical  of  the  radiation  heat  exchange  of  space 
craft.  Assume  that  the  surfaces  receive  the  solar  radiation  I,  and 
become  heated  and  in  turn  emit  long  wave-length  radiation,  which  is 
proportional  to  the  fourth  power  of  the  absolute  temperature.  The 
boundary  heat  transfer  is  then 


surface 


(51) 


where  a = Stefan -Boltzmann  constant 
e = surface  emittance 

If  the  inside  surface  of  a wall  of  finite  thickness  is  kept  at  a con- 
stant temperature  T , which  is  at  the  initial  temperature  when  t = 0, 
the  heat  transfer  relation  at  a time  t = t6  is 


Qt  = 


1 (Tt- 

i=0 


IT  ~ creV 


(52) 


T 

T 


must  be  found  from  the  following  relation  by  iteration 


Tt  Xq  + CTeTT4  - It  - (TT-i  " TQ)  (53) 

i=l 


Substituting  the  solution  T^  back  into  (39),  the  heat  flux  at  time 
t = t5  can  be  obtained. 
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Case  2:  Simultaneous  heat  and  mass  transfer  boundary: 

Transient  heat  transfer  through  a multilayer  solid  wall  when  one 
surface  is  wet  and  experiencing  either  evaporative  or  condensing  heat 
and  mass  transfer  is  treated  in  this  section.  The  surface  boundary 
condition  is: 

-x(|3)  = I + MTa  - T)  + KDL(Wa  - W)  (54) 

surface 


where 

I = Irradiated  heat  flux 

h = Convective  heat  transfer  coefficient 
c 

= Mass  transfer  coefficient 
T^  = Ambient  air  temperature 
= Ambient  air  humidity  ratio 
W = Humidity  ratio  of  saturated  air  at  the  surface 
temperature, T 

L = Latent  heat  of  evaporation  or  condensation 
Assuming  again  that  the  temperature  of  the  other  surface  of  the 
wall  is  kept  constant  at  the  initial  system  temperature  at  t = o,  the 
response  factor  relation  for  the  iterative  procedure  for  finding  T^ 
will  be  written  as 


Xo(tt-T0)  - ^‘c  Cra_'rT)  - KDL(Wa-WT)  = I - V X^T^-Tq) 

1 L-i 

i=0 
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As  can  be  seen  in  this  example , the  irradiation,  heat  and  mass  transfer 
coefficients,  air  temperature  and  air  humidity  ratio  can  also  be  treated 
as  time  variables. 

The  value  of  T must  be  found  first  from  (55)  by  iteration  and  put 
T 

back  into  the  heat  flux  equation  to  calculate  the  heat  transfer.  Many 

other  complex  heat  transfer  problems  can  be  solved  in  the  same  manner. 

10.  Calculation  of  Space  Temperature 

The  response  factors  developed  in  this  paper  are  not  only  useful 

for  evaluating  the  conduction  heat  transfer,  but  also  are  applicable 

to  the  calculation  of  a space  temperature. 

Consider  a simple  room  surrounded  by  a wall  whose  response  factors 

are  X.,  Y.  and  Z.  (j  = 0,  1, 

J J J 

Also  assume  that  the  room  air  temperature  at  time  t is  established 
as  a result  of  sensible  heat  balance  among  the  following  components: 


= heat  generated  at  time  T 
q^  = cooling  capacity  of  supply  air 
= 1.08  • (CFM)  • (T  - T ) 

T S , T 


where  (CFM)  = supply  air  flow  rate  in 
cu.  ft  per  min. 


T 

T 


= room  air  temperature  at  time 


T 


T 

S , T 


= supply  air  temperature  entering 
the  room 


q - heat  capacity  of  room  air 

d 


dT 

= V C = VC  (T  - T ,)/AT 

aadT  aaVT  t-1^ 
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where  V = room  air  volume 
a 


C = room  air  specific  heat 
a 


dT 

*i 

dT 


= time  change  of  room  air  temperature 


q = heat  gain  to  the  room  through  the  wall 
w 


= T . x.  + ) T . Y."1  A 
L Li  T-j  J L,  o,T-j  jj  w 


J=o 


J=o 


when  A = total  wall  area 
w 


It  is  assumed  that  the  response  factors  were 
previously  calculated  for  the  heat  flow  in 
the  direction  from  inside  to  outside  in  time 
increment  At 

The  heat  balance  equation  is 


qG  + qw 


~ qv  = 


The  room  air  temperature  T^  is  readily  obtained  by  substituting  each 

expression  of  heat  flow  into  the  above  equation  and  factoring  out 

T as  follows: 

T 


V c 
a a 

AT 


T 


T-l 


00 


CO 


[ Y T . x. 
L L>  t-j  j 

_i zl 


Y.  1 A + 1.08  (CFM)  T 

JJ  W s,T 


+ q 


G 


V C 

+ 1.08  (CFM)  + X A 
At  o w 
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The  knowledge  of  the  past  history  of  room  air  and  outdoor  air 
temperature,  therefore,  permits  the  evaluation  of  present  room  air 
temperature  with  the  use  of  response  factors.  Although  the  example 


cited  herein  is  a simple  one,  any  degree  of  complexity  can  be  added. 


if  necessary,  to  make  a complete  heat  balance  of  more  complicated 
system. 


According  to  equation  (29),  the  heat  flux  calculation  by  the 
convolution  relation  requires  a large  number  of  terms  for  the 


small  to  be  negligible.  When  the  heat  flux  values  are  to  be  cal- 
culated successively,  however,  it  is  possible  to  shorten  the 
computational  efforts  by  making  use  of  the  modified  response  factors. 
The  modified  response  factor  concept  may  be  explained  in  con- 
junction with  the  common  ratio  (CR)  relation  of  equation  (30)  as 
follows.  From  equation  (29)  the  heat  fluxes  at  r - r^  for  two 
consecutive  times  t-1  and  t can  be  expressed  as 


11,  Modified  Response  Factors 


summation  before  the  values  of  terms  ^ becomes  sufficiently 


CO 


(56) 


i=o 


00 


Ft  - = ) (X.  V.  . . • T Y.  V 
1, T~1  L r l , t -1 - l i 1 


i n,T-l-i 


) 


(57) 


i-o 
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By  multiplying  the  common  ratio  of  the  response  factors  (denoted  here 
“ B 6 

as  CR  = e P1  ) to  the  both  sides  of  equation  (57)  and  subtracting  it 
from  equation  (56), 


F,  " CR 

1 > T 


1,T-1 


<Xo 


V 


1,T 


r y 

o 


V 

n,T 


) 


00 


i=l 


CR)  V . - r (Y. 

1,T-1  1 


Y.  . . CR)  V [■ 
n,T-iJ 


By  noting  from  equation  (30)  that 


(58) 


X.  Y. 

= r—~  = CR  for  i > N + 1 


the  heat  flux  at  r = r^  for  time  t may  be  calculated  by 


N 


F1,t 

CR  • 

F1,t-1 

+ ) (x.1  v.  , - r y.'  v .) 

Li  l i,t-i  i n,T-r 

i=o 

(59) 

where  X . ' 

i 

= X.  - 

l 

Xi-1  • 

CR 

V 

= Y.  - 

l 

Yi-1  • 

CR 

(60) 

for 

i = 1, 

2,  ... 

N 

and  X ' - X 

o o 

Y ' = Y 
o o 
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• • o 


N),  are  called 


These  sets  of  finite  numbers,  X.',  and  Y.'  (i  = 0,  1 
the  modified  response  factors  of  the  first  kind. 

Since  the  value  of  N is  usually  around  15,  as  indicated  earlier 
for  most  of  the  building  walls  and  roofs,  the  calculation  effort  for 


factors.  The  similar  expression  can  be  derived  for  the  heat  flux 

at  r = r . , or  for  F . 

n+1  n,T 


Periodic  heat  flow 

The  response  factor  calculation  can  also  be  modified  to  shorten 
the  periodic  heat  flow  calculations.  Under  a periodic  boundary 
condition,  the  temperatures  and  heat  flux  must  assume  the  following 
relation 


where  p is  the  period  of  the  cycle  and  ft  = 0,  1,  2,  ...  <*>. 

Assuming  that  p is  larger  than  N,  beyond  which  the  response 
factors  are  evaluated  by  the  common  ratio  relation  (30),  equation 
(29)  can  be  expressed  as  follows: 


F. 


can  be  reduced  drastically  by  employing  the  modified  response 


for  o < r < p-1 
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whereby  X^',  Y^'  and  are  modified  response  factors  of  the  second  kind 
according  to  the  following  relationships 

X. 1  = X.  + X • (CR)i/(l  - CRP) 

x 1 p 

Y. '  = Y.  + Y (CR^/Cl  - CRP) 

1 x P 

Z. '  = Z±  + Zp  (CR)V(l  - CRP) 

12 • Cone lus ions 

General  formulae  for  calculating  thermal  response  factors  for 
multi-layer  structures  of  plane,  cylindrical  and  spherical  con- 
struction, have  been  developed  and  these  formulae  are  listed  in  this 
report.  Several  applications  of  these  response  factors  are  also 
illustrated  such  as 

a.  Interface  temperature  of  a multi-layer  construction 

b.  Evaluation  of  non-linear  boundary  temperature  problem 
such  as  radiation  and  evaporation 

c.  Evaluation  of  room  air  temperature 

The  computer  program  called  RESPTK  developed  to  obtain  the 
response  factors  based  upon  the  formulae  described  in  this  report 
is  found  in  the  Appendix. 


USCOMM-NBS-DC 
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Table  1 

Matrix  Elements  for  Plane  Layer 
Ay  = cosh  (qv£v) 

Bv  = K S(qviv) 

C = sinh  ) 

Dy  = cosh  (qv£v) 

where  jHv  = Vhl  ' *v 

K - 

* \) 


s(qv^v) 


sinh  (qvlv) 
v 


32 


m = 1 


Table  2 

Matrix  Elements  for  Cylindrical  Layer 


~ rv+l)  (^O  ,1  >2  + ^0*1^1  *2) 


- (<lv  \)+l^  (Ii  ji  Kq  + K1,1I0>a) 

where  I0jl  = IQ(qv^) 

lo , g ~ lo^Qv^+l) 

Ii  >1  — (iv 

52  = *1  (qv  *v+i) 

Ko,i  = Ko(qv^) 

*o,s  = K0(qv  ’v+l) 

Ki  5i  = Ki  (qvrv) 

^1  >2  = ^1  (qvrv+i> 


These  are  the  modified  Bessel  Functions. 
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Table  3 

Matrix  Elements  for  Spherical  Layer 


m 


\ ) (cosh(qv£v)  - -~T  S^qv^v)y 


rv  ' 


(rv+l 


bv  - ^ s(^v) 

a,  = 


JL  r 

'v  Rv 


! (qv2rvrv+l  - 1)  S(qvj?,v)  + cosh(qv£> 


Dv  = (-7^)  (cosh(qv£v)  + (— ) s(qv^v)y 


where  X y - r^+l  — 


S(qv^v) 


sinh(qvlv) 

qv^v 
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Table  4 

Derivatives  of  Matrix  Elements  for  Plane  Layer 

m = 0 


dAv  _ 

(bL'j  Sx  (q v£v) 

dp 

dBv  = 

l|“—  (qv^v> 

dp 

v2a/ 

H 

(^~  'j  TT"  r Si  (qvjfcv)  + cosh(qv£v)  j 
\2(V  Rv  L J 

dDy  = 
dp 

(*£) Si  (q*v) 

where 

Sj  (qv<v)  = Slnh<%^) 

qv^v 

q ^ \ = c°sh (qv^v)  - Sj  (qvxv) 

(qv^v)2 
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m = 1 


Table  5 

Derivatives  of  Matrix  Elements  for  Cylindrical  Layer 


}.)  (vIo,iKo, 


2 “ K0}1Io,2/ 


dBy 

dp 


= -(  rvr^±l )!  *1  »i  ^o  .a  ^3,  »1  Ip.a  \ + ( *V+.l 

' 2qv  K qvxv  J ^20v  /V  ^vXv 


dCv 

dp 


i-(-^^)(qvXv)(Io,lKl>S  + *0,1 


>2 


+ (-5±1  j(qvXv > i Kf) , 2 + Ki  »iIo,2  j 

2 C\> 


dDy> 

dp 


= (■£^±iXI°.lKo.=  • *0.1  ^.O  - - ^.,1..=) 


where  I0>1  , *o,3****^i>2  are  a H defined  previously  in  Table  2. 
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m = 2 


Table  6 

Derivatives  of  Matrix  Elements  for  Spherical  Layer 


dp  'ICAj' 


(I*ti  Sj(qv^v)  + (iy)  (fail)  s2(qvf,v) 


where  and  S2  have  been  defined  in  Table  4. 
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Table  7 

Limits  of  Derivative  Matrices 


When  (rv+i  - rv) /rv  is  sufficiently  small,  these  derivatives  can  be 


approximated  as  follows: 


dAv  = (iv3\  r,  i 

dp  '2a  ' L 2 \r v/  J 


dBv 

dp 


\2a  / L3 


+ (1-) 
12  Vrvy 


f_T\a 

lrvy 


dCv  = (l£\  f I 

dp  '2a  J L rv_!  Ry 


= (hi)  [l  + 3 
d \2a  / L 2 j 


38 


Table  7 (con't) 


dAv 

(h* 

) r2rv+i  + i) 

dp 

\6o^ 

/ \ rv  / 

dByj  = 

dp 

Rv  (■ 

Ty*i\(hL) 

rv  Aba»J 

dCy  _ 

L.  (: 

hi)  r2rv+i  + 2 /£vV 
2<V  L rv  3 \rv/ 

dp 

dDy  = 

(hL 

) (??)  fl?  * 0 

dp 

V60y 

When  (rv+i  - rv)/rv  is  sufficiently  small,  these  derivatives  can  be 
approximated  as  follows 


dAy 

dp 

. /v* 

2Qy 

(■  * i $ 

dBy 

dp 

lv3 

- 6<^ 

dCy 

= 

^v8T  £v 

dp 

Rv  1 

^ /L1  + rv 

dDy 

_ /.v2 

( i . 1 ivV. 

dp 

2aJ 

V1  + 3 rvA 
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m 


0 


B. 


v 


Table  8 

Matrix  Elements  for  p = for  Plane  Layer 

cos  E 

‘ Ry  sx 
E 


^ - 
dAv 

dp 

dBv 
dp 

dCv  _ 
dp 

dDy  = 

dp 

where 


Ry 

cos  E 


sin  E 


ysx  (E) 
yRySs  (E) 

1-  (S,  (E)  + cos  (E)) 
Ry 

■y  sx  (E) 

E 

Sx  (E) 


a. 

sin  E 

Sg (E)  - S^(E); cos--g 

E2 


y = 


„3 

J l V 

20y 
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Table  9 

Matrix  Elements  for  p = -g^  for  Cylindrical  Layer 
Av  - - f E2 (J0iY12  - Y01 J12) 

= | (-™Liy(J°lY°2  " Y01J02) 

Cv  ~ 2 Vrv+1y\"JllY12  + YhJ12;e| 

^ = ^ E2^j11y02  " Y1 1^02) 

dAy  _ TT_  f '| 

dp  4a  rrv+lry(JllY12  - Yll*Ji2)  + rv+l(JoiY02  " Y01J02>j 

d^~  " 4aVT^~y  {^rv+irv)(JllY02~YllJ02)  + (rv+l)S  (J01Y12’Y01J12)} 

df  = E2{^rv+lrv)(J01Y12-Y0iJi2>  + (ry+l)2  (J11Y02“Y11J02)'} 

dp  4a  {^rv+1  rv>  ("JQ1Y02  + y01j02)  " (rv+l)2  (~J11Y12  + Y11J12)| 
where  % = ^+1,  E1  = 7tf[  ^ 

J01  = -Jo  (Ei)»  J11  = Ji(Ex) 

J02  = «J0(E2)s  J12  - Jl  (Eg  ) 

K01  * Ko(Ei>.  KU  M (Ei  ) 

K02  = K0(Eg)s  K12  = Kx  (E2) 
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in  = 


Table  10 

Matrix  Elements  for  p - -3^  for  Spherical  Layer 


K = 


' V 


cos  E 


& V 

rv+l 


^ = YRv(5±IjSs(E) 

dp  x rv 

dCv  = y(lJ)(hL)a((*xE)dl  + i)si(E)  - (EjEg  + 1)S2  (E)) 

dp  Hr.Av  ^v2  ' 


= V(lv+1  Si  (E)  + (E)  j 


YV 

where  E = 

V 


v'  ' rv 


E,  = 


% rv,  E=  - []jk  S+l 


V TTT 


-v  T3^ 


Si  (E)  = jJ f - 


Ss(E) 


S,  (E)  - cos  E 


y = 25: 
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TO  (t) 


Sample  multi-layer  walls  of  Table  11 


Table  11 


layer 

V 

description 

& V 

(ft) 

(Btu/hr . ) 

Qv 

(ft2 /hr.) 

rv 

(ft) 

1 

Inside  air  film 
hx  = 1.20 

0 

— 

— 

5.000 

2 

Common  brick  4" 

0,333 

0.42 

0.019 

5 . 000 

3 

Face  brick  4" 

0.333 

0.77 

0.028 

5.333 

4 

Outside  air  film 
h0  = 3.0 

0 

— 

— 

5.666 
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Table  12 

Residues  at  p = 0 for 


0 < 

t ^ 6 

Dcp/B 

cp/B 

A Cp/B 

PW 

2.73117 

-2.94849 

7.88866 

CW 

2.6297 

-2.6874 

7.90340 

sw 

2.52953 

-2.44114 

7.89881 

6 < 

t ^ 26 

PW 

-2.31309 

3 . 36656 

-7.47058 

CW 

-2.19279 

3.07306 

-7.51776 

sw 

-2.07374 

2.79607 

-7.54387 

PW  = 

plane  wall 

m = 0 

CW  = 

cylindrical  wall 

m = 1 

sw  = 

spherical  wall 

m = 2 
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Table  13 

3ft : Roots  of  B(p)  = 0 


h 

PW 

CW 

SW 

1 

.17452 

.17701 

.17980 

2 

.84430 

.84634 

.848  6 6 

3 

2.56859 

2.57005 

2.57188 

4 

4.85967 

4.86146 

4.86360 

5 

8.85960 

8.86093 

8.86265 

6 

12.84988 

12.85127 

12.85303 

7 

19.15047 

19.15200 

19.15398 

8 

25.00846 

25.0095 

25.01083 

9 

33.33174 

33.33359 

33.33583 

10 

41.45064 

41.45137 

41.45249 

PW  = plane  wall 

CW  = cylindrical  wall 

SW  = spherical  wall 
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Table  14 

Response  Factors,  Btu  ft"  , F”  , hr~l 


i 

xi 

Yi 

Zi 

0 

PW 

.91949 

0.00013 

1.9834 

cw 

.92162 

0.00014 

1.97607 

sw 

.92363 

0.00011 

1.96864 

1 

PW 

-.16678 

0.00812 

-.51260 

cw 

-.16392 

0.00759 

-.52127 

sw 

-.16099 

0.00713 

-.52993 

2 

PW 

-0.07950 

0.03112 

-.23226 

cw 

-0.07744 

0.02916 

-.23749 

sw 

-0.07540 

0.02726 

-.24268 

3 

PW 

-0.05150 

0.04482 

-.15634 

cw 

-0.04987 

0.04185 

-.15997 

sw 

-0.04826 

0.03903 

-.16353 

4 

PW 

-0.03715 

0.04658 

-.11690 

cw 

-0.03580 

0.04340 

-.11954 

sw 

-0.03447 

0.04038 

-.12207 

5 

PW 

-0.02861 

0.04304 

-0.09216 

cw 

-0.02746 

0.04000 

-0.09410 

sw 

-0.02632 

0.03712 

-0.09592 

6 

PW 

-0.02292 

0.03784 

-0.07482 

cw 

-0.02192 

0.03508 

-0.07625 

sw 

-0.02094 

0.03247 

-0.07756 

7 

PW 

-0.01877 

0.03250 

-0.66173 

cw 

-0.01790 

0.03006 

-0.06277 

sw 

-0.01704 

0.02775 

-0.06369 

8 

PW 

-0.01556 

0.02761 

-0.05137 

cw 

-0.01480 

0.02548 

-0.05212 

sw 

-0.01404 

0.02344 

-0.05274 

9 

PW 

-0.01298 

0.02333 

-0.04294 

cw 

-0.01231 

0.02147 

-0.04346 

sw 

-0.01165 

0.01970 

-0.04386 

10 

PW 

-0.01086 

0.01965 

-0.03598 

cw 

-0.01028 

0.01804 

-0.03632 

sw 

-0.00970 

0.01651 

-0.03655 
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Table  14  (con't) 


i 

Xi 

Yi 

Zi 

11 

PW 

-0.00911 

0.01653 

-0.03018 

CW 

-0.00860 

0.01513 

-0.03039 

SW 

-0.00809 

0.01381 

-0.03050 

12 

PW 

-0.00764 

0.01389 

-0.02533 

CW 

-0.00719 

0.01269 

-0.02544 

SW 

-0.00675 

0. 01155 

-0.02547 

13 

PW 

-0.00642 

0.01167 

-0.02126 

CW 

-0.00602 

0.01063 

-0.02131 

SW 

-0.00564 

0.00965 

-0.02127 

14 

PW 

-0.00539 

0.00980 

-0.01786 

CW 

-0.00505 

0.00891 

-0.01785 

SW 

-0.00471 

0.00807 

-0.01777 

CR 

PW 

0.8398 

CW 

0.8378 

SW 

0.8358 

PW  = plane  wall  m = 0 

CW  = cylindrical  wall  m = 1 
SW  = spherical  wall  m = 2 
CR  = common  ratio  where  for  i > 15 


Xj+1  _ Y i+l  ^ i+1 
xi  ~ 


CR 
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t 

24 

23 

22 

21 

20 

19 

18 

17 

16 

15 

14 

13 

12 

11 

10 

9 

8 

7 

6 

5 

4 

3 

2 

1 


Table  15 

Plane  Wall  Model 


TIt 

T0t 

Qlt 

c$t 

Exact 

Solution 

Response  Factor 
Solution 

Exact 
Solut ion 

Response  Factor 
Solution 

75 

77 

-17.15 

-17.15 

31.04 

30.58 

75 

79 

-19.20 

-19.18 

32.92 

33.77 

75 

81 

-20.90 

-20.83 

39.47 

38.71 

75 

83 

-21.48 

-21.33 

44.79 

47.35 

75 

85 

-20.00 

-19.91 

69.48 

70.30 

75 

87 

-17.10 

-17.10 

-20.15 

-24.28 

75 

138 

-13.72 

-13.76 

-72.88 

-78.01 

75 

162 

-10.45 

-10.53 

-101.46 

-103.27 

75 

168 

-7.82 

-7.91 

-111.45 

-111.14 

75 

163 

-6.06 

-6.13 

-98.60 

-99.37 

75 

148 

-5.05 

-5.07 

-76.34 

-74.90 

75 

128 

-4.04 

-4.36 

-32.97 

-32.06 

75 

104 

-3.86 

-3.89 

-28.51 

-28.87 

75 

100 

-3.65 

-3.68 

-23.32 

-23.04 

75 

95 

-3.73 

-3.75 

-16.17 

-16.18 

75 

90 

-4.13 

-4.15 

-10.75 

-10.70 

75 

86 

-4.86 

-4.87 

-4.76 

-4.73 

75 

82 

-5.82 

-5.83 

7.48 

8.04 

75 

76 

-6.94 

-6.95 

14.46 

14.34 

75 

74 

-8.21 

-8.22 

17.11 

17.42 

75 

74 

-9.65 

-9.66 

19.08 

18.77 

75 

75 

-11.29 

-11.30 

19.94 

20.27 

75 

76 

-13.11 

-13.12 

25.02 

24.84 

75 

76 

-15.07 

-15.08 

27.90 

28.39 

TI:  Inside  temperature,  (F) 

T0:  Outside  temperature,  (F) 

QI:  Inside  heat  flux,  (Btu  hr" ^ f t”^) 

Q0:  Outside  heat  flux,  (Btu  hr  ^ ft”2) 


CO  CO 


= 4YTIt-j 

j = 0 

j ~° 

CD 

CD 

Q0  = ;Y.-TI  . 
*t  J t-j 

- )Z. -T0 
— J 

j = o 

j = o 
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Cylindrical  Wall  Table  16 


t 

Tit 

T0t 

Qit 

Q0t 

Exact: 

Solution 

Response  Factor 
Solution 

Exact 

Solution 

Response  Factor 
Solution 

24 

75 

77 

-17.94 

-17.95 

31.68 

31.23 

23 

75 

79 

-20.13 

-20.12 

33.76 

34.62 

22 

75 

81 

-22.00 

-21.88 

40.37 

39  .80 

21 

75 

83 

-22.60 

-22.44 

46.17 

48  .76 

20 

75 

85 

-21.06 

-20  .96 

71.40 

72.21 

19 

75 

87 

-18.01 

-18.01 

-17.72 

-21  .91 

18 

75 

138 

-14.44 

-14.48 

-75.48 

-75.61 

17 

75 

162 

-10.99 

-11.07 

-99.28 

-101.12 

16 

75 

168 

-8.21 

-8.30 

-109.69 

-109  .39 

15 

75 

163 

-6.33 

-6.41 

-97.30 

-98  .08 

14 

75 

148 

-5.25 

-5.29 

-75.50 

-74.04 

13 

75 

128 

-4.50 

-4.53 

-32.40 

-31.48 

12 

75 

104 

-3.99 

-4.02 

-28.06 

-28.43 

11 

75 

100 

-3.76 

-3.78 

-23.02 

-22.73 

10 

75 

95 

-3.82 

-3.85 

-15.97 

-15.98 

9 

75 

90 

-4.24 

-4.26 

-10.64 

-10.59 

8 

75 

86 

-5.00 

-5.01 

-4.73 

-4.70 

7 

75 

82 

-6.00 

-6.01 

7.47 

8.04 

6 

75 

76 

-7.18 

-7.18 

14.49 

14.36 

5 

75 

74 

-8.50 

-8.51 

17.18 

17.49 

4 

75 

74 

-10.02 

-10.03 

19.23 

18  .9  2 

3 

75 

75 

-11.75 

-11.75 

20.17 

20.50 

2 

75 

76 

-13.67 

-13.68 

25.36 

25.18 

1 

75 

76 

-15.75 

-15.76 

28.37 

28  .8  6 

TI: 

Inside  temperature,  (F) 

T0: 

Outside  temperature,  (F) 

QI: 

Inside  heat  flux  (Btu  hr--*" 

ft"2) 

00: 

Outside  heat  flux  (Btu  hr--*- 

ff2 

Wt  - h Tit.j  - r IYj  T0t-j 


QO  = >Y.  TI  . - T )z.  T0 

u J t-j  L J ^t-j 


5.666 

5.000 
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Spherical  Wall  Table  17 


Time 

TI 

w~ 

QI 

Q0 

(hr) 

Exact 

Solution 

Response  Factor 
Solution 

Exact 

Solution 

Response  Factor 
Solution 

24 

75 

77 

-18.74 

-18.75 

32.21 

31.80 

23 

75 

79 

-21.07 

-21.05 

34.52 

35.39 

22 

75 

81 

-23.03 

-22.95 

41.60 

40.82 

21 

75 

83 

-23.74 

-23.57 

47.49 

50.12 

20 

75 

85 

-22.15 

-22.04 

73.26 

74.06 

19 

75 

87 

-18.94 

-18.93 

-15.32 

-19.57 

18 

75 

138 

-15.18 

.-15.22 

-73.10 

-73.24 

17 

75 

162 

-11.54 

-11.62 

-97.13 

-98.98 

16 

75 

168 

-8.60 

-8.69 

-107.96 

-107.65 

15 

75 

163 

-6.60 

-6.69 

-96.03 

-96.81 

14 

75 

148 

-5.46 

-5.50 

-74.68 

-73.20 

13 

75 

128 

-4.66 

-4.69 

-31.86 

-30.92 

12 

75 

104 

-4.11 

-4.13 

-27.64 

-28.02 

11 

75 

100 

-3.85 

-3.88 

-22.74 

-22.45 

10 

75 

95 

-3.91 

-3.93 

-15.80 

-15.82 

9 

75 

90 

-4.32 

-4.35 

-10.57 

-10.52 

8 

75 

86 

-5.11 

-5.12 

-4.74 

-4.71 

7 

75 

82 

-6.16 

-6.16 

7.42 

7.99 

6 

75 

76 

-7.39 

-7.39 

14.46 

14.33 

5 

75 

74 

-8.78 

-8.78 

17.20 

17.51 

4 

75 

74 

-10.37 

-10.38 

19.32 

19.01 

3 

75 

75 

-12.18 

-12.19 

20.33 

20.67 

2 

75 

76 

-14.21 

-14.22 

25.62 

25.45 

1 

75 

76 

-16.41 

-16.42 

28.77 

29.27 

TI:  Inside  temperature,  (F) 

T0 : Outside  temperature,  (F) 

QI:  Inside  heat  flux,  (Btu  hr”^  ft"2) 

Q0:  Outside  heat  flux,  (Btu  hr-^  ft”^) 


00  CO 


Qit  =ZvTlt-j 

- ryY..T,  . 

J t-j 

t_i. 

n 

o 

o 

II 

00 

00 

v*1 

- )Z,.T . 
J t-j 

j = o 

r _ / 5 . 666  ,2 
\5 . 000/ 

o 

II 
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m = 1 . 


Table  18 

Formulas  for  Cylinder  and  Sphere 


G'  = 


dG'  = 
dp 


m = 2 


dp 


= Ai  \ 

I? 

Jx  (Ex  ) 

vri J 

J0(Ei) 

'£i! 

.TiiSL 

\rj\ 

/Lej0(e) 

1 

cos  (E)\ 

In  A1 

L - 

S:  (E)  J 

'rx2 

. r cos 

U A 

^2a 

/L1  ‘ 

at  p = i0 


(E)  - 


J0(E) Jx (E) 


n 


Jo^  (E) 


Sx2(E) 


where  definitions  of  E,  Sx  (E)  and  S2 (E)  are  identical  to  those  used  in  the 

previous  tables.  As  Ex  approaches  to  zero  G'  approaches  zero  for  both 

2 \ 


1 a o dG'  . \(\i\  , . , /ri  \Ai 

” - 1 and  2>  if  becomes  2a  for  m = 1 and  -I  — A~ 


—J  for  m = 2. 
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Table  19 


Response  Factors  for  Semi-infinite  Region 
Common  Symbols 


L _ n and  n ^ 6 


-@2  i/fJL 


O 02 

0 ^ = (—)  I — - — - — — dB 

eJ{Yo2(0)  + Jo2 (e)3 


m = 0 Z,  = 2 L ' ^ 

1 V tt 

\ = \ (/2  - 2) 

z.  = z,  (A  - 2 /in  + /m> , i > 3 

1 

m = 1 Zx  = I+x  (0X ) 


Z2  — IjjL  (02  “ 20^  ) 

zi  = Lp,  (0 j_  - 20i_1  + 0i_2)  for  i > 3 


m = 2 

Z1  = 2L 

/t 

V tt 

(1  + 

At 

2 /i j,/ 

\ = 2L 

/H” 

V tt 

(/2 

- 2) 

Zi  = 2L 

nr 

'V  tt 

(A 

- 2/m  + /T-2)  for  i > 3 

These  relationships  show  a very  interesting  fact  such  that  Z^ 
i = 2,  3,...®  are  identical  for  the  cases  where  m = 0,  1 and 
from  the  cases  of  m = 0 and  1,  it  should  follow  that 

(01  - 20!.!  + 0t_2)  = ^j=-  (fi  . 2/TT  + fPZ) 

2\2  f°  1 - e'^1^  dB 

a)  J0  e3uc2  <e>  + j02  (P)J 

which  seems  to  be  a remarkable  relationship. 


for 

2.  Moreover, 
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APPENDIX 


Response  factor  formulas  developed  in  the  main  text  of  this  paper 
were  used  in  the  computer  program  called  RESPTK.  The  Fortran  listings 
of  RESPTK  and  other  necessary  subroutines  to  calculate  the  thermal 
response  factors  of  various  multi-layer  heat  conduction  systems  such 
as  depicted  in  Figure  4 are  attached  herewith.  The  main  program  to 
perform  the  input/output  operation  for  RESPTK  is  called  RESP.  Sample 
input  and  output  for  RESP  obtained  for  the  systems  described  in  Figure  4 
are  also  attached.  Certain  portions  of  the  computer  program  are  written 
in  Fortran  V (Univac  1108)  and  certain  modifications  to  the  program 
will  be  necessary  for  use  with  a compiler  that  does  not  recognize 
statements  made  for  the  Univac  Fortran  V compiler. 
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Input  Requirement  of  the  Computer  Program  RESPTK 
RES PIK  (K,  L,  R,  G,  AG,  KG,  X,  Y,  Z,  NL,  DT,  NR,  CR,  U,  IM,  IS,  F) 

Input : 

K = Thermal  conductivity  (BTU/hr,  ft,  F)  of  each 

layer  given  in  the  order  for  minimum  radius  to 
the  larger  radii  (Fig.  2).  For  the  plane  wall, 
it  should  be  given  from  inside  surface  layer  to 
the  outer  layers.  For  the  layer  with  no  thermal 
mass,  such  as  surface  boundary  layer,  conductance 
values  should  be  used. 

L = Thickness  of  each  layer  (ft)  given  in  the  order 
for  minimum  radius  to  the  larger  radii  (Fig.  2). 
This  could  be  zero  for  some  layers,  i.e.  surface 
boundary  layer . 

R = Radius  (ft)  of  each  layer  boundary  given  in  the 
order  of  minimum  value  to  the  larger  values 
(Fig.  2).  For  plane  wall  model,  any  arbitrary 
value  being  same  for  all  the  layer,  should  be 

given.  Note  that  the  number  of  R is  NL  + 1. 

2 

G = Thermal  diffusivity  (ft  /hr)  of  each  layer  given 
in  the  order  for  the  minimum  radius  to  the  larger 
radii  (Fig.  2).  For  the  layers  with  no  thermal 
mass,  such  as  surface  boundary  layers  and  air 
space,  G should  be  zero  for  this  program. 
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AC 


KG 


X,  Y,  Z 


NL 


DT 


NR 


Thermal  diffusivity  (ft  /hr)  of  soLid  core 
or  semi- inf inite  layer  (Fig.  1),  given  only 
when  IS  = 1,  or  2. 

Thermal  conductivity  (Btu/hr,  ft,  °F)  of 

solid  core  or  of  semi-infinite  region 

(Fig.  1),  given  only  when  IS  = 1 or  2. 

2 

Thermal  response  factors  (Btu/hr,  ft  ) 
generated  by  this  program  for  a wall  of 
finite  thickness. 

Number  of  layers  to  be  considered  for  the 
heat  conduction  system.  Surface  boundary 
layers  and  air  spaces  should  be  treated 
as  separate  layers  with  G = 0 (Fig.  2). 

When  IS  = 1 or  2 , the  solid  core  or  the 
semi-infinite  region  should  not  be  counted 
as  a layer . 

Time  increament  (hr)  for  which  the  calcu- 
lation of  heat  flux  is  desired.  For  hourly 
calculation  DT  = 1. 

Number  of  X,  Y,  and  Z generated  by  the 
program.  NR  is  the  output  of  this  program 
such  that  the  values  of  X,  Y,  Z can  be 
calculated  by  a common  ratio  CR  as  follows: 


X(J  + 1)  Y ( J + 1)  = Z(J  + 1) 

X(J)  Y(J)  Z(J) 


when  J ^ NR. 
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CR 

U 


IM 


IS 


F 


Common  ratio  described  above. 

Overall  thermal  conductance  obtained  by  the 
reciprocal  of  total  thermal  resistance  of  the 
heat  conduction  system  under  consideration, 
Btu/hr,  ft^,  °F. 

Curvature  index  (Fig.  1) 
if  IM  = 0 plane  system 

- 1 cylindrical  system 

= 2 spherical  system 

Heat  conduction  system  index  (Fig.  1) 
if  IS  = 0 finite  wall 

- 1 semi-infinite  region  attached 

= 2 solid  core  attached 

Response  factors  for  the  system  with  solid 
core  or  semi-infinite  region. 


Calculation  of  heat  flux 

(A)  Referring  to  Fig.  2,  the  heat  flux  QI(N)  and  QO(N)  can  be 
evaluated  as  follows,  where  N is  the  time  index  such  that 
time  = DT*N„ 

M M 

r-*  -—i 

QI(N ) = ^ X(J)*TI(N-J)  - GM*  ) Y(J)*TO(N-J) 

J = 1 J = 1 

M M 

QO(N)  =)  Y (J )*TI(N-J ) -V  Z(J )*TO(N-J ) 

J = 1 J = 1 
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2 

where  QI  and  QO,  and  TX  and  TO  are  heat  fluxes  (Btu/hr,  ft  ) 


and  temperatures  (F)  at  surfaces  where  the  radii  are  minimum 
and  maximum,  or  at  the  inside  and  outside  surfaces.  The  values 
of  QI  and  QO  are  positive  when  heat  is  flowing  from  TI  side  to 
TO  side  or  from  inside  to  outside. 

In  above  equation  for  QI(N), 

M = maximum  number  of  response  factors  to  be  used,  value  of 
which  will  be  determined  by  the  significance  of 
X(M)*TI(N-M).  Usually  M < 72  (when  M > NR,  X(J),  Y(J) 


and  Z(J)  should  be  calculated  by  the  common  ratio  CR 
such  as  described  earlier)  and 


_ R(KL  + 1) 
R(L) 


**  IM,  which  is  unity  for  the  plane 


wall  problem. 

B.  For  calculating  heat  conduction  for  the  system  with  the  semi- 
infinite region  (when  IS  = 1), 


M 

QI(N  ) = > F(J)*TI(N-J) 

4-J 

J = 1 


TI(N-J)  is  the  temperature  at  the  surface  where  the  radius  is 
minimum  (inner  surface)  at  time  (N-J)*DT.  The  value  of  QI(N) 
is  positive  when  heat  is  flowing  in  the  direction  from  the 
minimum  radius  (inside  surface)  to  the  larger  radii  (to  outer 
layer  and  to  the  semi-infinite  region). 
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(C)  For  calculating  the  heat  conduction  for  solid  core  system 


(IS  = Z), 

M 

QQ(N ) = -)  F(J ) TO(N-J)  , Btu/hr,  ft2 

Lj 

J = I 

TO(N-J)  is  the  temperature  of  the  surface  where  the  radius 
is  maximum  (outside  surface).  The  heat  flux  QO(N)  is  positive 
as  it  is  defined  in  the  above  equation  when  heat  is  flowing 
in  the  direction  from  the  smaller  radius  to  the  larger  radii. 

Bessel  function 

The  calculation  for  IM  = 1 (cylindrical  system)  requires  a 
double  precision  Bessel  function  subroutine  in  the  following 
forms: 

JQ(X)  = DBEJ (X,  0) 

JX(X)  = DBEJ (X,  1) 

Yq(X)  = DBEY(X,  0) 

YX(X)  = DBEY(X,  1) 
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ZU0  KEAD(b»2)  L ( 1 ) * K l I ) f Ij  ( X ) * C ( I ) » REb  ( 1 ) 

1 F U N • eCJ  • z:  • AND  • I M • EO  • U ) oO  TO  .601 
READ  K t RHO  t AND  U OF  GROUND  IF  lN=l 
BOO  IF ( In.NE.o)  Re AD l S » 2 ) Kb » UG » CG 

Ab  ThERMal  DIFFUSIVITY  OF  EaRTH 
IF ( In. NE • 0 ) Ab— Kb/Cb/Ub 
IF (NLAYR. Eb.O)  Go  TO  Sol 


i f ( i i‘t . t-0 • u ) uu  ro  jo  j 
Kt.Alj  lbf?i  ( 1 ) » i - 1 » NUL  A Y h ) 
ou  1 u bOL 
OUl  K(l)=lu. 

Uu  bOb  I — ^ f (jNLh  Y r< 

h ( 1 ) ~K  l I - 1 J + L ( A ) 

0 i i 2 i F ( 1 1 j » fc_  W . ^ • MUU  • I i*i » i4L  • U ) k£  aL)  ( b * 1 1 2 J ( K ViK  b l J ) » J-l  » 4 ) 

Ub  1 1 3 I — x » Nl . A Y R 

1 1 3 ( K I ■! t\  ( if  J ) > J — 1 » 4.  1 

it-  ( Iixj.lO.  x ) KtAlMbfli^)  (RW.Kbtu)  » J-l»4) 

LU  1 U9  Il1»NLaYR 
A <-  ( L ( I ) ) li(t  » ill  f lilt 
All  b ( 1 ) — tj  » 

h { I ) ~ 1 t/nLbl  1 ) 

0 u T u 104 

Ail)  b ( I ) =K  ( I J / L ( 1 ) / L)  ( 1 ) 
a 1 1 9 LuNT  anUF 

bU  i b|vjA~  l H l NuLAYR  ) /K  l 1 ) ) **iM 
V1;  K 1 1 t ( t>  f c 0 ( ) 

£<j7  FuRF'aI  12mA  ) 

b mLL  K tbi-  I K ( K »L»h  »u»au»Ku»X»Y  */*uLaYN»U'LlTaT  » I iK  T » CR  * U 1 » 1 M » I h » F ) 
WKlTt,(b»lUU) 

IF  ( 1 i*i  • t W • u ) m k a T t_  l b » 7 o A ) 

/ U 1 FoWMATlSur.U  RLAufc  WALL  bYSTF.M 

1 F J I )*i  « LO  • i J WK  1 Tt  ( b r ) 

/U2  FoKMmT  Ibunlt  C Y L 1 UUR i C aL  bYSTtM 
IF  ( I lb  • LO  • a ) wK  i T i_  ( b f 7 Uo  ) 

I U3  FoRFiMTlbUnU  ShHrkiCML  bYbTFM 
A h I T r ( b » 1 U i ) 

ViK  1 f r.  ( ti  * auc  ) 
h h I T t ( h m bu ) 

|M|  i l_  Y K ♦ LU  . I)  ) v.iO  T U hll^: 

it-  ( 1 N • e.O  • L • A UU  • I iv>  • i Ah.  • 0 ) «*  K i T F ( U » A 2 U ) KbMJb*Lto»  l R ibK  b ( U ) f J—  i * 4 ) 

L'U  2U?  I — a » N L A Y K 

it-  i L l i ) ) ^u2  f 2bb  » fUt? 

L 1 1 3 K ( 1 ) — ()  . 

lu2  VtKlTh  lb»  ibO)  1 * L ^ I ) * i\  l i ) » I ) l I ) » L ( i ) » KUb  ( i ) t ( RMK  ( I ► J ) » J=1  » 4 ) 
i F ( 1 1 ) . t.  9 . 1 ) m/KiTc.  ln»  L 2 u ) Kb*l)0»Lb»  iKMKbi  J)  » J—  1 t 4 ) 
a 2 0 h uKNiaT  l It- ^7  . 3»  iF'iO  .2  t lh  10 . b»  10X  » 4 Ad) 
bu2  Wk  I Th.  ( b » I Ub  ) UlL  I A 1 
AhITfcl(b»  a u 4 ) U 1 
to  K I T t£  ( b » i u b ) 

Ink  i ic.  ( b » 4U U ) 

IF  ( lU.IMt  . 0 ) bO  To  ibbb 
A K I T L ( b * A U 7 ) 

UO  1 14  N-A  t NWT 
OIm  = N-  1 

1 14  WKlTblrsflutl)  JN  * X ( N ) f Y I N ) » Z l N ) 
bo  Tu  r>()4 
lubb  ink  IT  c.  ( o * nbb  ) 

Dbb  t- uhI',iaT  ( buiiu  U F 

it-  ( 1 i\t  • lO  • i ) bb  To  bUb 
it-  ( Il>j  .tr.o  . L > AivlO  • IM  • to  . 0 ) bO  TO  bOb 
LU  hub  NLir|\jr<7 
U|\j  — N—  1 
X ( N ) L-a  ( Im  ) 

DUb  wRlTc.(o»bOb)  JI>i»X(i4) 
bO  Tu  bOn 

D U b Uu  b04  (\j—  i t NRT 
U t j — I g — 1 

b09  Ah  i Tt  l b f nub  ) Jim»(-In) 
oU8  h uRibM  T ( 1 a f 1 h 2 1 • b ) 


:3  J 4 w k I T ( h # 4 <J  li  ) 

Vv  K i T i . ( t > * 4 U U ) 
wk  I Th.  (n*  x L / ) Ck 
U (NTLST.r.fci.O  ) 
CALL  ) it. A I t.  a t i f / » * 


TO  .Ail) 

I * To  » uLL  1 mT  # NH  » nK  I » CM  A » Lk’  ) 


bo  I U 0 0 U 

• 

oOO  STOP 

LNI) 

■ : 


. 


r,  r 


b I 1 FOR  b Ai-iPlE  » bAMFL  L 

7 Hi  S PROGRAM  .1  LLUbTK  A 1 t-b  I HE  USE  uF  RLSFTk  AND  HEATX 
Ue  VElOPE.u  uY  7 hE  NaTIOnAl  oURtAU  Oh  STanuAKUS 
C I i iNbiuE  7 lMPEK  A 1 UkE  WHERE  K lb  T ML  bMALLEbT 

l fu  Ou  I b 1 1 )h  7 lMPEK A ( UkL  WHERE  R ib  1 ML  LARGEST 

0 R 1 Ht-KhiAl.  L OhDUC  T i v I I Y , R ib  THE  THRMAL  CONDUCT  AhCF  OF  1 HE 

0 L A Y h (\  IF  i r ifc  kE  lb  NO  I IiEkKmL  MaSS  bUCIt  A b AlK  SPACE 

C L IHiLKNE.bS  oF  ThE  LAYER  L~(J  IF  NO  THERMAL  MASS 

C F kAl/IOS  OF  i HE  LAYER  USE  aRdITRaKY  VALUE  FOR  (HE  PLANE  WAl-L 

L b IHtKHAL  UlFFUSl  ViT  i Of-  7 ME  LmYEK  6—0  FUR  NO  THERMAL  MASS 

u i ME  Mb  I Oi  j N ( 1 u ) , 1_  ( .1  0 ) , K ( 1 0 ) , (■>  ( 1 0 ) , X ( 1 0 U ) , Y (10  0)  f Z ( 1 0 0 ) 

LiMEnSIOm  t.V(AUtlj  ,b0L(20G)  » T L T ( 2 0 0 ) ,FU0u)  ,00(200)  ,TZ1200)  ,TU(200) 

L>  i Mt.  mS  i On  iXI'+o)  » T Y ( a m ) , hTvv  ( 1 1 ) 

o i MEub ION  ii  (2  0 0 ) » i 2 ( 2 U 0 ) ,02(20  01 


2 F OK M A T l 1 u x / ) 
b h uR  M a 7 ( 1 u h /«h) 

i u4  h oRMa  T ( hLihU  kiiALL  COMHObl  I 1 UN  SIaIEo  F RUM  INSIDE  SURFACE 
xUb  FuRMaT ( l2uh0  LAYER  R L 

1 rr  L 

iiO  FORMAT  thijnb  PLaNl  imAi_L  MODEL 
j-il  FoRMaJ  IbUr.O  C Y L l nDR  I C ML  vnALL  MODEL 
i i 2 F 0Rf“t  a7  (Sumo  SFHEfIlAl  wAlL  MODEL 


) 

) 

) 

) 

) 


L 

0 

o 

c 

c 


Kc  Al  R » L 

REaD(S»2i  n D h Y , i I M A X 
i F ( Nu  A Y • t_u  • 0 ) OO  To  10U 
1/  ROOM  »uK  I L -VERA  i UkE 
UN  UUTuUOK  AiR  T EMRc.n A I UkE. 

o Ni>u  VE.LOlITY  In  rNoTS 

S o L b o l A K kAuIaTi On 

7j7  TlIaL  LlOuD  AmuuM 

Kt.AD(b,l)  (T/(o7  1,J7  — J»24) 

Kfc.  AD  ( b » 1 ) ( Dm  ( OT  ) , UT-  1 , 2*4  ) 

R L A D l S , 1 ) ( w V ( J I ) * J I - I » 2 4 ) 

RE  Ai ) ( b * 1 M bOL ( jT  ) , oT= i » 24 ) 

RE  All  ( b » 1 > ( I C T ( oT  ) , lT-1 , 24  ) 
l EoRMHT(l2ho«u) 

DO  1 Ob  II=l,m0aY 
OU  iub  jK-1^24 

0 f =24* l I I ”1 ) 4 JR 
Urs  l J I )— Dtt(>jK) 

Il;(  Jl  ) -Drt  ( uR  ) 

W V ( J I ) -W  w ( uK  ) 

ToT(JT)=lOl (JR) 
boL ( uT ) =boL (JR ) 
i U3  I c ( J I ) = TZ ( oR ) 

JIUO  READ  (b,  2)  ImL»IM 

K Cl  A D ( S , 5 ) ( K ( I ) , i = 1 , N L ) 

RhAD(S,5)  ( L ( I ) , i = 1 » NL ) 

nll=nl+ 1 

K E A D ( b » b ) (R(l)»i=l»NLL) 
KeAD(S,S)  ( o ( 1 ) * 1 = 1, Nl) 
RtAI)(b»b)  uT 
j.  h ( 1 M • E.  Q • O ) wRiTE(b,llU) 

1 F t 1 ivi  . t.0  . i ) WR 1 ft  ( 6 , 1 i i ) 

IF  ( i M • E u • 2 ) W K I "I  t ( 6 * 1 1 2 ) 

Wk1TE(h»104) 

wrITcl  (o,  iub) 

DO  4 J-l , ULL 

4 toklTe.(o’4>  K(J),l(J)»R(J)»6(D) 
■S  h ORMA T ( 6F2U  • k ) 


A o=0 
Rb=0 


iN=0 

l/X-H  l NLL  ) /k  ( 1 ) 

it-  ( I H . t.  W • U ) i . 

1 F ( 1 i"’l  ,t.^t  It  OM  A — u X 
1 F ( I M • L W • '■■■  i It  I'!  A — <.T  X ♦ ♦ 8 

C/\LL  UlSL  lK(KrL»rt»b»rt^*Ko»A*  Y » / * hL  * U V » NK  » CK  » U » 1M  * I U » F ) 

H | W ( <i  ) — 0 • / 
n i mi  ( 4 ) =0  • y 
h | W ( ui  ) - i . 

HTW ( 1 ) =D I 
hi  I W ( ) — l\)K 
h I W ( S ) -48 
hi T W ( 4 ) — <vi /-a 
h [ W ( b ) ^Cr< 
to  K 1 1 fc_  ( n t i . 8 ) 

13  huHMATU^unl  DHlKlJ  7/MhT)  SOL(KT)  TCT(KT)  to 

1V(*T)  FOe  10  1 MtATtoT  ) 

018=0  . 

0 8 H = U . 

UO  11  KT-49 » NTwAA 
C aLL  F u ( ;v  \>  i K r ) » L » F OC  » F u T » 1 ) 

ti  | w (tv)  —Foe 

h i W ( / ) -SuL i K T ) 

01  jh  — |J(  \ ( K t ) 

h i w ( 1 1 > — 1 1 i (KT  h 

IN  fHlb  HKUuKAM  lHt  PRlSLN i I I Nth  lb  t>l*48  iH  HUUri 
L/U.L  ktV  » ( I 0»  1 X »4H#K  1 )r 

o a Ll  i<  c.  V I l I L • F f » h H f K i ) 

0 mLL  HL  A I A h X » Y t Z.  • TX  * i Y » OuB  * T oT  » Ht_  A I to  I » Hi  to  » o IF* » 02P  ) 

While.  ( 8 » i tL  ) l Jh  h K I ) # T / I h 1 ) * bOL  (t\T)  » I C 7 C K T ) » to  V ( K I ) » F UC  » 1 O T » Ht.  a I W T 

I 8 F o H M a T l 1 o F 1 2.  • 8 ) 

1 1 ( K I -48 ) = I 0 I 

I 8 ( K I -48  ) — I 8 ( K 1 ) 

O A ( K I -48  ) -HLA'f  to  T 

II  To(K  i UTu 

toh  I l L ( / ) I i. 

to  K I T l l 7 ) 1 L 
toj<  I Tc.  ( 7 ) 

KF  to  1 Hi , 7 
Li  jh 


i>)l\i  F OR  A » A 

buBRoU  I Inc.  RESk  f K ! K » l » k AG  » Kb  * X » Y »Z»  NL  » D1  » NR  » CR  » U » 1 M » I S » F ) 
DIMENSION  k ( 1 U ) » l ( 1 U ) * k ( 1 0 ) » b ( 1 0 ) » A ( 1 Ou  ) » Y ( 100  ) » Z ( 1 00 ) t AP  ( 10  ) » tiP  ( 1 
10)»CH(10)»uP(lu)»A(lU)»bUU)»C(lU)rUllU)»ZRl(i)  » ZR2  ( 3 ) »RB(3)  »RAP(3 
2 ) » RUOT ( iuu ) ( 3» IU0 ) »ZRk (b? 10U ) >RX ( 10U ) »RY ( 100 ) * AZ ( 100 ) 


3 i K ( 1 U U ) 

REAL  K»L»KG 
P I ~3 « In  1 

1 F ( I b » E.O  • 1 • aND  • i M • E Q • 0 ) to  R i,  T (r  lb»oO*+) 
oU4  F OR mm T ( SonU  bEMl  — XnF  iNITE  PLmNE  WALL  ) 

it-  ( I S • E O » i • a NU  « I M • E Q • 1 ) rtRiTE(fc>»oOb) 
oiiS  F o R M h T l 5 u r i u bEM l™ INFINITE  CYLINDER  ) 

IF  ( I b . lQ  . i • AND  ® I m * EQ  • 2 ) *R1 T£ (b»  t>0o  ) 

o06  FuR Mat (sumu  bEMi- infinite  bPHEKE  ) 

IF  ( I b o EG  • l • And  • I iW  ® EO  • 0 ) tfiR  i Tt  ( b » t>0  7 ) 
otl7  FoRMa T(5uhb  b O L I L)  SLaH  ) 

i F ( I b s t.0  • • AND  •Ii'#teEQ*it  WkITE(6»6U2) 

oU2  F ORMaT  ( 5brt  bOL  I U-- — C 1 L 1 NITER  ) 

IF  ( Ib»t.Q«ir,»ANL)*  I M • E Q • 2 ) wRiTE  (b*  o(Jb ) 
ou3  FoRMaT  (5uFi  SoL  I D-”--"”""’  bPFtERc.  ) 


M3=3 

i F t I b a tW  e • AND  » X m o NE  • 0 ) 1*1  3 = 1 

i F ( I b « NE * i ) bO  To  b!3 
OUH  £L=Ro/W ( NL+1 ) 

UY=R  t N L -*■  1 ) **2/aG/UT 
CaLL  faRFlUt f ZL» lM»AZ  ) 

JLF  d b • t.0  • 1 ® AND  • Nl  • EQ  « U ) b O TO  Pul 
oi3  CaLL  AbCu2 ( 0 • * R t L * R r b » aX » HX * C X » DX » XM * NL ) 

Kh ( 1 ) — U X 
Rh<2)=l . 

KH ( 3 ) — AX 
U = 1 . /bX 
uo  1 1=1 t nL 
PX  = U 

CaLL  AbCuP2 (PX  »K ( I ) » L ( i ) » R ( I ) t b ( i ) * AP ( 1 ) » RP ( I ) » CP ( I ) » DP < I ) » IM ) 
1 OmLL  aBlU*(PX»MI)  >L(i)  >R(D  »b(l)  r a ( I ) »o(I)  »C(I)  »D(D  rlM»l) 

1 F ( Nl . L T « l i ^0  To  b02 

OhLL  DERV I (A*b*C*Q» AP»bP»CP»l)P»APP»bPP»CPP»DPP»NL) 
ou  To  bOo 
su2  ApP=Mp(l t 
bRp=f^p  ( l ) 

CpP=oP ( 1 ) 

DPP=DP ( 1 ) 

bU3  lF(Ib.N£.2)  bO  To  bOl 
iF  ( I Fi » E 0 • O J bO  To  bOl 
CaLL  boL  iu  i 0 • » R ( 1 ) »Ko»AGpIM»  HF  » F(FP  i 
ZKl  ( I ) = (-lkP+FIFP*AX)  /DX/oT 
LK2( 1 )=-^Kl ( 1 ) 

WKlTc(brlbuO) 

lbUO  FORMAT  1 12UFI  CPP  HFP 

1 aX  DX 

1 *+  U 0 FoRMIaT  (4F^0.S) 

WK  ITt ( b» 1 AL0 ) CPP»HFP* AX»DX 
bO  TO  12 1L 
bUl  KAPU)=DkP 
RmP(2)=0. 

RAP ( b) =AFF 

UO  2 1 = 1 » b 

C 1 =R  AP ( I )/bX/DT 

C'd  = Rb(  i ) *bPP/BX/HX/DT 

ZR2 ( 1 )=»C14C2 

’ 2 ZRl ( i )=-ZR^( I )4Rb( I 5 /BX 


1 e:  1 2 W(<1  Tt  ( b*o*+  ) 

b4  l-  uRMm  r 1 5)uriu  RLSiuuES  mT  P=0 
WKHtthriDU)  ( ZR  1 ( 1 ) » I = X r ivi  3 ) 

VVHlTt(tr>piUO)  ( 2R2  ( I ) , I — 1 t M,}  ) 
iuO  P oRMa  T (3P2U«b) 

C ROOTS  of-  HlP)=u. 
c.\?.  IMN,  A X — 4 U 

II-  ( I b . LG  . 2 . A iXD  . I !*i . nL  . 0 ) nMaX=1  ill) 

Px  = 0 .001 
URU  — u . 1 / u I 

It-  <Ib.tW«2.  AimD  . I N • iMh  . 0 ) uPO=3.141b*3.l4lo*Mb/|<(l>/R(l)*0.2b 
bi.X^tj . UOu  1 

if-  < Ib.tO.c  .AND.  I M • HL  . 0 ) l)L_X— l)PO/  1 UOU 
N=o 

A A I T t ( b 1 t>  O ) 

b3  P oRNa T l bUnU  koO  i 5 0|-  Li  ( P ) — 0 

11  Ul_=l'PU 

CmLL  MhCU/r'  l P x r t\  ri_  t K t (j  * pX  r 8X  » C X * L)A  * lM  * NL  ) 

it-  ( 1 S • tQ  • 2 • AND  • In  • NL  • 0 ) u ALL  SuLuX  ( PX  * k ( 1 ) »Kb*AG»IM»HX*DX»TE5Tl) 
lb  PxP  — P X +()L 

LaLl.  AtiCUt.  t P xp  » K » L t R » b * A XP  » bXP * C XP  * L)XP  » I lvl » NL  ) 

1 p ( I "> . nE  . 2 ) i-jO  T u 2 1 0 
if-(lM.tO.U)  bO  TO  (r'Li 

LmLL  bOLUA  ( P X P » K ( 1 ) »t\b.A o * iM»HXP»OXP»  TclS  12 ) 
il-  (TbSI  1*  lt_ST2)  Il2*il3»il4 
1 1 4 P X — P xP 

T LST  l =T  Eb  1 
bo  To  lb 

a 12  IP  ( 0 l — uLa)  130»lOU*ll7 
x 1 7 UL.=UL/2. 

oo  To  i b 

113  lf-(TtSll)  ilH»lly*llrt 
HP  KxX  = PX 

bo  To  bl 
11H  K X X — P X P 
bO  TO  Ol 

130  AH  = At5S  l Ttb  I 1 /TLb  I 2 ) 

KXX=(PX+MD*PXP) / ( 1+An) 
bu  To  31 

213  lp(HX*bXpt  12  » 1 3 » 14 
14  pX=PXP 
bX=HxP 
bO  To  lb 

12  Ip  ( Oi_-OLx  ) oO  * 30  * i 7 

17  UL=DL/2. 
bO  TO  15 

13  i F ( H x ) 1 « > 1 9 * 1 8 

19  HXX=PX 

bo  TO  31 

18  kxX=PXP 
bO  To  31 

30  AB=AHS (Bx/bXP ) 

KXX=(Pa  + AD*PXP)  / i 1 • + A ti  ) 

31  N=N+1 

HOOT  (Im)  =KXX 

if<n.gt»d  dpo=woot ( n ) -root (n-i ) 

nrt=n 

WpITL(b*4l)  N*kOOT(N) 

41  p oRMaT l 1 10  * lF2u . b ) 

PX=R  x X +Di_A 
1 EST  MX— 4o 
TpSTx=KXx*oT 


IF  (TFS1  X-'ftSTMX>4?r42»43 

42  IF  (N.Ll  .nMmx  ) GO  TO  11 

43  WR  I Tt  ( t> » o5  ) 

r>S  F oRMAT  ( 5lih0  RLSUDUES  mT  P=ROOT(N)  ) 

L/ b«J0  Jo=jLH\|RT 
P*=ROO") ( Oo ) 

DO  Si  u=l *nL 

L aLL  ^yCUtlPXfKt  j)  * L ( J ) » K ( J ) » G ( J ) » a ( J ) » B ( J ) » C ( J ) » 0 ( J ) » I M » 1 ) 
hi  Call  abCuP2 (px » k * j) *l i j) »r ( j) »g( j) » ap (o) >bp( J) »CP( J) »op ( j) » im) 

CaLL  AhCu2  IPX»!\»l*R»o»aX»hX»CX»OX»IM»  NL  ) 

IF  (NL.LT.  oO  To  bi)4 

CaLL  DtRw'l  (A»ti»C»D»  AP  * t>P  »CP»  DP  » APP  * BPP » CPP  » DPP  t NL  ) 

Go  TO  50  b 
hl)4  ApP— ah l 1 ) 
bPP— np l 1 ) 

CpP-LP ( 1 ; 

OpP=UP ( 1 ) 

ou5  ip(Ib.NE.z)  *40  To  214 
Ip(lM.tO.O)  GO  TO  214 

CaLL  SoLiulPX»K(i)  » Kb  » AG  » I h » HF  » HF  P ) 
if-  ( HR  ) 4ui  » 4l)0  » 4dl 

auI  PYS  = InF *aX-Ca ) /Pa/PX/ ( UPP-HFP*bX-HF *bPP ) /UT 
Ou  To  4 02 
4 U 1)  PyS=u. 

4 0 2 Rm(1»JJ)=r'IS 
GO  To  fiO  1 

4 1 4 PI =BPP*P**PA*U I 
Km ( 1 f OJ ) =uX/HY 
Ka<2# JJ)=1./PY 
K A ( 3 » Jo ) — AA/PY 
biJl  p/=Px*uT 
32  kx ( JJ ) =Eah ( -p2 ) 

5 K Y ( Jo ) - ( 1 . -ExP (. Pc ) ) **2 
ollO  WRlTr.(b»b4)  ROOT  ( JJ  ) . ( RA  ( M » OJ  ) t M=1  r M3  ) 
b4  FoRMaT  14I-2U  .b  ) 

DO  1 b4  Ou— 1»NR1 
OO  lh4  M = i » M 3 

ZrI (M)=Rm lM» JJ) *RX ( Jj) + ZK1 (M) 

1 b4  zK2  ( M)  =Ra  l i*l  r JJ)  *kX  ( Ju)  *RX  ( JJ)  * ( 1-2/RX  (oJ)  ) + ZR2  ( M ) 

11=1 

111=2 


80  FoRMaT (buHU  RESPONSE  FACTORS  OF  FINITE  SLAB  ) 

bl  FORMAT l l2Uh0  J X(J)  Y(J) 

1 Z(J)  ) 

/ui  format(  12ufi1  response  factors  for  solid  cylindrical  objects 

i ) 

/02  FoRMaT (12Uhl  RESPONSE  FaCTORS  FOR  SOLID  SPHERICAL  OBJECTS 

1 ) 


IF  ( IS.  Eft  • 2 . AND  • I M • EO  el)  wRlTE(br701) 
1 F ( IS.LQ.2.A  NO*  I Ri.EQ.2)  uRlTE(b>702) 
IF ( IS.dO.UJ  wR1TlTO»8U) 

Wr I TE ( b » toi ) 

IF<Zk1 (2; .lT.O)  Z R 1 ( 2 ) = 0 • 

WRlTt(orbb)  II» (ZRl(M) »M=1»M3) 
hRlTE(b»b3)  I 1 1 r ( ZR2 ( M ) * M=  1 » M3  ) 

DO  b 7 M= 1 * M3 
ZRK ( m » 1 ) =ZKl ( M ) 
b7  ZRK (M»2)=ZK2(M) 

55  FORMAT ( 1 1 0 » 3F2U . b ) 

N T = 1 U 0 

DO  58  N=3»Nl 
NR=N 


UO  hi  M=i»!*i3 
o 1 Z K R ( Ri  » N ) = u • 

UO  6 1 M-  i * M3 
UO  57  oO=i»NkT 
PV  = (KX (GO) ) **N 

b7  ZkK(i4»N>=ZRK(M#N)+PZ*HY ( J J ) *R A ( M » JJ ) 

KK I Tfc ( n * bb ) IM  » l ZKK  l M » N ) t 16=1  * M3  ) 
ip (N .L I . h)  Go  10  58 
1 1-  ST1=/Rk  U » ivi ) / Zr<K  ( 1 *iM-l  i 
Tr  5T2=ZRk  1 1 * N-i ) /ZkK  l i * N-2  ) 
lpSTb=ARb(  lE6Tl-tEbT2) 

IF ( TESY  3-u . OOOul ) b9*b9»5R 

b8  CONTINUE 
ti9  UO  hi)  I'J—  i t NK 
X ( N ) -/KK lirN) 

1 (N)-/KMc:)N) 
hU  Z ( N ) = / K K Uj  * N ) 

LK=TES1 2 
whITt.  ln»h^)  Ck 

h2  PoKMmTIIuHU  Ck=1F10.6) 

ih  (iS.hO.t  • At4U.  l6.EQ.0i  OO  TO  HUU 
it-  ( 1 n.NE.  i ) GO  TO  OOU 
yijl  ip(NL.hQ.G)  bU  To  005 

bi-^^Kb/bwKT l U I *rtG*P i ) 
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UO  204  JLbUMvlR 

tU4  a/7g)=oF* (SORT (Zo)-2.*bOKT(ZJ-l . ) +SuRT ( / G-2 • > ) 

N K R — < •J  R 
CO  Tv)  30  u 

olU  uu  3ti  1 J-NKrhO 
Z ( J+i )-/(^)*OR 
X(0+J  )=X  lc)*CR 
3U 1 Y ( J+ 1 ) = Y ( U ) *cR 
NkR — 60 

jlJ 0 UO  2_l)b  J-i  » NhR 

cU6  p ( J ) -X  l J t - Y (u)*Y(G)/(ZlG)  + AZ  l J ) J 
t — Nkk 
CO  To  9 On 
905  UO  9U4  G=i>NR 
9U4  P ( J) -A/ lo) 

9u6WKlTt.(h»/o7) 

lu7  K)KMaT  (buhl)  G 

C K 1 — 1 * 

UO  20 8 G=1 » 50 

CK=P ( J+l ) / E (G) 

1 p_STCR  = Abb  (CK-CRi  ) 

IP ( TEST Ck-O • OOUO 1 ) 611»611*612 

ol2  Ckl=CR 
o J— G— 1 

ZU8  WKlTt.(h»ZU9)  JO  ft-  ( G ) 

Zt)9  FORMAT  ( liiU»  1F20.5) 
oil  NR =0 

Cr=Ck 1 
bo  TO  90 U 
oGO  wkITpJ(to»2u7) 

UO  210  0=1 t NR 
p ( j ) =2#Y (o)-(X(J)+Z(o) ) 

OO— 0—1 

zlO  WRlTt(hrZU9)  OO  M-  ( O ) 

9U  o return 
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1 ) » C T l 1 u ) * f / 1 l I 0 ) » a f T ( i o ) » b T I ( ] 0 ) » C T I ( 1 u ) # U T I (10) 
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o ! ( o ) -oP  v o t 
U I < o J -UP ( o J 
2 oOl'iT  iNUt 
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L>oUpu_p_  PkllISIoN  LMbu » ubt Y » Zw 1 » Z02 

Kt ML  K ( 10 ) »L ( 1U ) » JOl » ou2» Jll * J12 
P i = 3 . l4lbV27 

i *0  . b 

ip  ( NL  • LT  • <£.)  K ( 2 ) — K ( 1 ) + i_  ( i ) 

UQ  4 1—1 »uL 

ip  (b(  I ) ) iti3»  1 U 3 * 1 1)2 

x U2  ip  (2  ) i t i » 101 
x U 1 2lv=SwR  I ( /_/  O ( i ) ) 
i-  U 1 — 2 Q * R i i ) 

Zu 2-/_<  j*R  i j.  + i ) 

Zw  L — ZO  + L ( x ) 
i P ( I >i « i\iP  . x ) DO  To  /> 
o 0 1 — . ) i i c.  J l »-A>  1 » U ) 
oil  = JhtJ  i 1 »1) 

01)2— UBP..J  l *_u.2  » U > 

0 i 2 — i)ptJ  l 4o2  • 1 ) 

Y U 1 — I )r(t  Y Ix(j1»U) 

1 ll=obt_Y  WiOl » 1 ) 

Y u2— i iht  Y 1 x_o2  » u ) 

V 1 2 — i j ) i t.  Y l O 2 * 1 ) 

AX(  I )--PP*^2*  ( Jul*Y  12- Yu  1*0 12) 

bx  < 1 >-PP*w  U + l ) /t\(  1 ) ♦ (-Y0 1*002+ Jul*YU2) 

Lx ( I ) =K ( x > 2K  l 1 + 1 ) * (-0  1 1*  I 12+Yll*ol2)  *Pp-*ZQ2*/02 
Da ( 1 ) =PP*2u2* (*jll*Y02  — Y ii*oU2) 

Ou  To  4 

3 Lu-S  x l\i  l 2oi_  ) 

L l=LuS  ( 2wi_  ) 
j1-Cu//Cvl 

b2= (bl-Cl ) /ZoL/ZuL 
i p ( I Pi  • t 0 • ) bU  T u b 

AX ( I ) =L  1 

bX ( 1 i=L ( i )/K ( i ) *bl 
LX<  I )=-7ul*K( 1 )/l( i ) *uu 
Da ( 1 ) -L 1 
bU  T u *4 

b obi— K (1+1)2 K(i) 

MX  ( 1 ) = oto*  l Ll-L  (I)/K(i+i)*Sl) 
bx  < 1 )=•_(  i ) /K  ( I ) *ol**Sl 

Lx  ( 1 )=L(  i ) *L(  i ) /K  ( i ) 2K  ( I ) *r\(  I )/L  l I )*(-(Zol *Zo2  + 1 ) *S1+C1 > 

bO  To  4 
l AX ( 1 ) — 1 • 

C X ( I ) — U • 

UX(I)  = IRU  + 1)/K(a)  )**1M 
iP(iPi.tW.u)  rfX(l)=L(  i)XK(I) 
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GO  TO  4 

103  A x ( I ) = 1 • 

L5  X ( 1 ) = 1 / IN  ( 1 ) 

Cx  < I ) =u . 

UX(I)  = (R(i  + l)/K(i>>**ilvl 

4 LONT i NUt 
a-AX  ( 1 ) 
b— BX ( 1 ) 

L=CX ( 1 ) 

U=DX l 1 ) 

ip  (Ni_.LT  .2  ) oO  TO  b 
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toll  l tOR  U t I ) 

SviMKOU  f I Nr.  AhCuPtf  ( Z » K » L * K » 0 » AP  # HP  t CP  » Df  » lM  ) 
bulIhLL  Pr<r_GlSioN  ZG 1 » Zu2  t Dfcjt  J » Urtt  Y 
KtAL  K»L»Uul»JU2*Jll»Jl2 
K 1=.} . 141 
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/Li,  1 — Zu  *K 
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ir  ( I ,Ji . oh  « j ) oC  To  .■> 

X — H * ( R + L J 
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/.  1 = ( r<+L  ) /r 

Util  — ilntJ  It.ul  » l.l  ) 
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011  — Ui  it  J l /-G  1 » i ) 

012  — Out.  J l 4-  G 2 # 1 ) 

Till—  l )tit  Y 1 ^ G 1 » il  1 

V l)2  — urit.Y  1 ZG2  » 0 ) 
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t i2— i jrtt Y l z g2  » 1 J 
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la'-  ( a ♦ U2+f(t  1*j()«l)-Y*(-J11*Y.L2+Y11*J12)  ) *PP 

GO  To  h 
A x-;L  *(_♦  U • :■>/  G 
K j.  — K +L 
Kr.S— L/r\ 

Co— S i N ( Zgl. ) 

C l = CuS IZgC ) 

Sl=Cu/ZOL 

1 1 1 I i*i  • t Q • b 1 i-»0  To  b 

Ap  = **(Kl*bi/K-L*'jP/KI 
bP-RtS*X  *Ki  *S2/R 
ml =/ol 
Zr2=/ajc 

Ct  = X*  (i_/r<  ) »*2/Kt;b*  ( (2.»K»ki/L/L+l  ) *S1-(ZP1*ZP2+1.  ) *S2  ) 

(kl/K*bl  + li_/K)*(Rl/K)  *S2 ) 

GO  To  4 
b Ar=X*S1 

CP=X  * (bl+41  ) /Ht.S 

U f — X * S i 
GO  T J 4 
1U.5  /R=U. 
b H — t ) . 

CR  = 0 . 

L)P=U  . 

GU  To  4 

1U1  it  ( l.'I.NF  • << ) GO  To  b 
X=L*L*U . b/G 
At  — X 

bP  = X +L/K/ b 
CP=K/L*X*^# 

UK  — X 
GO  To  4 

b it  (lb.  NE  • 1 ) oO  TO  7 


H 1 -H  +-L 

AH=  ( U * .1  ) i ( i<  1 /H  t ) *u  • 5/G 
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bR  = ti . h/G* i u • >->* l Ki*R 1-k*R ) 4k l/R-R*Rl*LOo ( K l/R ) ) 
IjU  i v>  4 
7 a— l.*L*u  • ."•»/' o 
h i — K +L 

A h - X / 3 • * t ci  4 R | / t<  f J • ) 

|)p-L/K»ki/h*x/j. 

U-'-fv/L^X+L/H+L/'^f^.+lUKl/L/l  +ll  ♦ nbooh?  > 
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4 hr  1 1 JKiM 
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I'UllbLf  lJKLC  islu  J I JmF>J  » /_bo 
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b U T U 3 0 u 

a 0 0 C-COb ( ZUi ) 

b-Si.iizbi ; /zbl 
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Pl2=2./Pi 
th=O.UUl 
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wkITt(n»xuu) 
tok  I TE ( o » lu A ) 

xul  kOKMATlSonU  K Z(K) 

Z i 1 ) =2*Zi_*bUkT  llj)  /bQTPi 
Ll-l l 1 ) 

Z l 2 ) -Z ( 1) * ( buk I ( z . ) -z . ) 

GO  Z K-^OU 
Zl\  — K 

2 Z IK  ) ~Z  ( 1 ) * IboRl  ( ZK  ) -z .*buRT (7K“1 > +bQRT  (ZK-2.  ) ) 

IF  ( To  VtJ 

1 F <1 N . t G • 1 ) Tu  1 

L.  I 1 ) -Z  ( 1 ) + ZL 

GO  T o 71) 

1 A — P I Z *LUO  ( I)  • b*cO  ) *0 . 3o74bF>9l 
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1 A = G 
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bu  1 / L - 1 > b 0 » < U 
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Z jO  = unL  J l zG  r i) ) 
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1 a - 1 u 
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GO  28  L=1»lT 
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28  Zb(L)=(l.-zP)*ZTlL) 
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T I MF  INCRFMFNT  DT= 

l. 

THFRMAL 

CONDUCT  ANCF 

11=  1.002 

RFSPONSF  FACTORS 

j 

0 

F 

.73488 

i 

? 

— •POO  38 

-.  toons 

3 

4 

-.07084 

-.05563 

h 

b 

-.04593 

-.03869 

7 

-.03284 

6 

-.02797 

9 

-.02385 

10 

-.02035 

u 

-.01736 

IP 

-.01481 

13 

-.01264 

14 

-.01079 

COMMON  RATIO  CR-  .85332 


SPHFhlCAL  WAI.I 


L A f I - f<  L (I.  ) K ( I ) (T)  C ( T ) RFS(T)  DESCRIPTION 

No  OF  LAYFRS 


1 

.0(H) 

.000 

.00 

.000 

.92 

IMS TOE  SURFACE ( 

2 

. o ) n 

. 130 

Ao.no 

.30  0 

. no 

1 / A— IN  1 TNOl  FUM 

3 

.500 

1.000 

i a n . n n 

.200 

.00 

G-TM  CONCRETE 

.700 

120.00 

.230 

GROUND 

TIME  INCREMENT  DT= 

l . 

THF  RMAL 

CONDUCTANCE 

U=  .514 

. R FSPONSF  FACTORS 


J 

E 

0 

• 890  45 

1 

-. OAOftl 

2 

-.02555 

3 

-.01379 

4 

-.00841 

5 

-.00533 

fS 

-.00344 

7 

-.00223 

ft 

-.00145 

9 

-.00093 

10 

-. 00050 

l 1 

-.00038 

12 

-.00024 

13 

-.00015 

14 

-.00009 

15 

-.00005 

1 <3 

-.00003 

17 

-.00002 

1ft 

-.0000] 

19 

-.0i)001 

20 

-.00000 

21 

-.00000 

22 

-.00000 

23 

-.00000 

24 

-.00000 

25 

-.00000 

«?.n 

-.00000 

27 

- . 0 0 000 

?,H 

-.00000 

29 

-.00000 

COMMON 

RATTO  CR-  .80254 

SPHERICAL 


LAYER 

NO 

L ( I ) 

K ( I ) 

( I ) 

C ( T ) 

RES ( I ) DESCRIPTION 

OE  LAYFPS 

.250 

50.00 

.320 

8-TN  RAn  PAPFR 

l 

.uin 

.130 

fto.on 

• 300 

.00  1/8- IN  LINOLEUM 

2 

• non 

.000 

• on 

.000 

.92  SURFACF 

TIME  INCREMENT  DT= 

1 . 

THERMAL 

CONDUCTANCE 

U=  1.002 

R F SPO  N SC  F ACTORS 


J 

F 

0 

.71441 

1 

— • 2 3 0 4 A 

2 

-.12057 

3 

-.08440 

4 

— • 05333 

5 

-.04880 

ft 

-.03758 

7 

-.02909 

A 

-.02254 

9 

-.01747 

10 

-.01355 

1 1 

-.01050 

12 

-.00814 

COMMON  RATIO  CR=  .77522 


